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

Amin Ebrahimi  Department of Materials Science and Engineering, Delft
University of Technology, Mekelweg 2, 2628CD Delft, The Netherlands  Chris R.
Kleijn  Department of Chemical Engineering, Delft University of Technology,
van der Maasweg 9, 2629HZ Delft, The Netherlands  Ian M. Richardson

###### Abstract

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éclet number.

##  1 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 ($\vec{S}_{\mathrm{d}}$) is
added to the momentum equation,

$\vec{S}_{\mathrm{d}}=\frac{-\mu}{K}\vec{V},$ (1)

where $\vec{V}$ is the fluid velocity vector and $\mu$ 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_{\mathrm{L}}$) as follows:

$K=\frac{\mu}{C}\frac{f_{\mathrm{L}}^{3}}{\left(1-f_{\mathrm{L}}\right)^{2}},$
(2)

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}\,$\text{\,}\mathrm{kg}\text{\,}{\mathrm{s}}^{-1}\text{\,}{\mathrm{m}}^{-3}$$
[21, 22]; however, values between $10^{4}$ and
$10^{8}\,$\text{\,}\mathrm{kg}\text{\,}{\mathrm{s}}^{-1}\text{\,}{\mathrm{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 $\mathcal{O}(10^{4})$ [23] and $\mathcal{O}(10^{8})$ [24] have been
proposed for stainless steel, and values between $\mathcal{O}(10^{4})$ [25]
and $\mathcal{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 phase-change 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 (CaCl2H12O6) 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 enthalpy-porosity 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.

##  2 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 ($\mathrm{D}=0.1~{}$\mathrm{m}$$). The
enclosure is initially filled with solid material at a temperature
($T_{\mathrm{i}}$) below the melting temperature $T_{\mathrm{m}}$ (for
isothermal phase-change) or solidus temperature $T_{\mathrm{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_{\mathrm{i}}$) to the hot wall temperature ($T_{\mathrm{h}}$),
whereas the cold wall is kept at a temperature $T_{\mathrm{c}}$ =
$T_{\mathrm{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 $\mu$. 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_{\mathrm{h}}$, and the right vertical wall is kept at an isothermal temperature $T_{\mathrm{c}}$, whereas the horizontal walls are adiabatic. Table 1: Thermophysical properties of the artificial materials used in the simulations. Property | Value | Unit
---|---|---
Density $\rho$ | $10^{3}$ | $\mathrm{kg}\text{\,}{\mathrm{m}}^{-3}$
Specific heat capacity $c_{\mathrm{p}}$ | $10^{3}$ | $\mathrm{J}\text{\,}{\mathrm{kg}}^{-1}\text{\,}{\mathrm{K}}^{-1}$
Thermal conductivity $k$ | 1, 10 and 100 | $\mathrm{W}\text{\,}{\mathrm{m}}^{-1}\text{\,}\mathrm{K}$
Dynamic viscosity $\mu$ | $10^{-3}$ | $\mathrm{kg}\text{\,}{\mathrm{m}}^{-1}\text{\,}{\mathrm{s}}^{-1}$
Latent heat of fusion $L_{\mathrm{m}}$ | $10^{5}$ | $\mathrm{J}\text{\,}{\mathrm{kg}}^{-1}$
Thermal expansion coefficient $\beta$ | $10^{-6}$ | ${\mathrm{K}}^{-1}$
Melting temperature $T_{\mathrm{m}}=\frac{\left(T_{\mathrm{l}}+T_{\mathrm{s}}\right)}{2}$ | $1.5\text{\times}{10}^{3}\text{\,}$ | $\mathrm{K}$
Melting-temperature range $\Delta T_{\mathrm{m}}$ | 0, 5, 10, 25 and 50 | $\mathrm{K}$

##  3 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 $\mathrm{Fo}=t\alpha/\mathrm{D}^{2}$,
$\vec{V}^{*}=\vec{V}\,\mathrm{D}/\alpha$,
$T^{*}=T/\left(T_{\mathrm{h}}-T_{\mathrm{melt}}\right)$,
$H^{*}=T^{*}+f_{\mathrm{L}}/\mathrm{Ste}$,
$p^{*}=p\,\mathrm{D}^{2}/\left(\rho\alpha^{2}\right)$ and
$K^{*}=K/\mathrm{D}^{2}$ [10], the unsteady momentum and energy equations are
cast in conservative dimensionless form, respectively, as follows:

$\frac{1}{\mathrm{Pr}}\left(\frac{\partial\vec{V}^{*}}{\partial\mathrm{Fo}}+\vec{V}^{*}\cdot\nabla\vec{V}^{*}+\nabla
p^{*}\right)=\nabla^{2}\vec{V}^{*}+\mathrm{Ra}T^{*}\hat{y}-\frac{\vec{V}^{*}}{K^{*}},$
(3) $\frac{\partial H^{*}}{\partial\mathrm{Fo}}+\vec{V}^{*}\cdot\nabla
H^{*}=\nabla^{2}T^{*}+\frac{1}{\mathrm{Ste}}\left(\frac{\partial
f_{\mathrm{L}}}{\partial\mathrm{Fo}}+\nabla\vec{V}^{*}\right),$ (4)

where the Prandtl number (Pr) represents the ratio of momentum diffusivity to
thermal diffusivity $\left(\alpha=k/\left(\rho c_{\mathrm{p}}\right)\right)$
in molten materials and is defined as

$\mathrm{Pr}=\frac{\mu}{\rho\alpha}.$ (5)

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

$\mathrm{Gr}=\frac{g\beta\rho^{2}\,\mathrm{D}^{3}\left(T_{\mathrm{h}}-T_{\mathrm{melt}}\right)}{\mu^{2}},$
(6)
$\mathrm{Ra}=\mathrm{Gr}\cdot\mathrm{Pr}=\frac{g\beta\rho\,\mathrm{D}^{3}\left(T_{\mathrm{h}}-T_{\mathrm{melt}}\right)}{\mu\alpha}.$
(7)

The Stefan number (Ste) is the ratio of sensible to latent heat and is defined
as

$\mathrm{Ste}=\frac{c_{\mathrm{p}}\left(T_{\mathrm{h}}-T_{\mathrm{melt}}\right)}{L_{\mathrm{m}}}.$
(8)

In the above, $\rho$ is density, $\mu$ dynamic viscosity, $p$ the static
pressure, $k$ thermal conductivity, $c_{\mathrm{p}}$ specific heat capacity at
constant pressure, $f_{\mathrm{L}}$ the local liquid fraction,
$L_{\mathrm{m}}$ the latent heat of melting or solidification, $\beta$ the
thermal expansion coefficient, $g$ the gravitational acceleration, $t$ time
and $\vec{V}$ the fluid velocity vector. $\hat{y}$ is the unit vector in the
y-axis direction. $T_{\mathrm{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:

$f_{\mathrm{L}}\left(T\right)=\frac{T-T_{\mathrm{s}}}{T_{\mathrm{l}}-T_{\mathrm{s}}};T_{\mathrm{s}}\leq
T\leq T_{\mathrm{l}},$ (9)

where $T_{\mathrm{l}}$ and $T_{\mathrm{s}}$ are liquidus and solidus
temperatures, respectively. In the case of isothermal phase-transformation, 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.,
$\nabla\vec{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].

$\frac{-\vec{V}^{*}}{K^{*}}=-\frac{\mathrm{D}^{2}\,\vec{V}^{*}}{K}=-\frac{C}{\mu}\frac{\left(1-f_{\mathrm{L}}\right)^{2}}{f_{\mathrm{L}}^{3}+\epsilon}\,\mathrm{D}^{2}\,\vec{V}^{*},$
(10)

where $C$ is the permeability coefficient and $\epsilon$ 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 $\left(f_{\mathrm{L}}=1\right)$ while its
limiting value for $f_{\mathrm{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 $\theta$ of melting-temperature range
$\left(\Delta T_{\mathrm{m}}=T_{\mathrm{l}}-T_{\mathrm{s}}\right)$ to the
temperature difference between the hot and cold wall $\left(\Delta
T_{\mathrm{w}}=T_{\mathrm{h}}-T_{\mathrm{c}}\right)$ is defined as

$\theta=\frac{T_{\mathrm{l}}-T_{\mathrm{s}}}{T_{\mathrm{h}}-T_{\mathrm{c}}}=\frac{\Delta
T_{\mathrm{m}}}{\Delta T_{\mathrm{w}}}.$ (11)

The volumetric fraction of liquid inside the computational domain at a time
instant $t$ is evaluated as

$\phi(t)=\frac{1}{\mathrm{W}\cdot\mathrm{D}}\iint\limits_{\mathrm{domain}}f_{\mathrm{L}}(t)\mathop{}\\!\mathrm{d}A.$
(12)

We use the relative difference
$\Delta\phi(t)=\left|\phi_{2}(t)-\phi_{1}(t)\right|/\phi_{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
($\Gamma$) is employed to quantify the difference between the solid–liquid
interface morphology predicted using different permeability coefficients
$C_{1}$ and $C_{2}$ as follows:

$\Gamma=\frac{1}{\mathrm{W}\cdot\mathrm{D}}\int\limits_{y=0}^{y=D}\sqrt{\left(x_{\mathrm{f},C_{1}}-x_{\mathrm{f},C_{2}}\right)^{2}}\mathop{}\\!\mathrm{d}y,$
(13)

where $x_{\mathrm{f}}$ is the position of the liquidus line.

##  4 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
$\Delta x_{\mathrm{i}}/\mathrm{D}=\Delta
y_{\mathrm{i}}/\mathrm{D}=$4\text{\times}{10}^{-3}\text{\,}$$ 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., $\left(\Delta T_{\mathrm{m}}/\Delta
T_{\mathrm{w}}\right)\mathrm{W}\leq\Delta x_{\mathrm{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
$\left(\Delta x_{\mathrm{N}}=\Delta
x_{\mathrm{i}}/2^{\mathrm{N}},\mathrm{N}=1,\cdots,4\right)$ were applied every
ten time-steps. A fixed time-step size of $\Delta\mathrm{Fo}=k\Delta
t/\left(\rho c_{\mathrm{p}}\Delta x_{\mathrm{i}}^{2}\right)=10^{-4}$ was
selected, corresponding to a Courant number $\left(\mathrm{Co}=|\vec{V}|\Delta
t/\Delta x\right)$ less than 0.25.

The conservation equations were discretised using the finite-volume method.
The central- differencing 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 second-order 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 $\phi$ from one
iteration to the next is less than $10^{-10}$.

##  5 Results

###  5.1 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\text{\,}\mathrm{K}$ above the liquidus temperature $T_{\mathrm{l}}$ =
$919\text{\,}\mathrm{K}$. The phase-change process starts at $t$ =
$0\text{\,}\mathrm{s}$ by suddenly changing the temperature of one side of the
slab to $573\text{\,}\mathrm{K}$, below the solidus $T_{\mathrm{s}}$ =
$821\text{\,}\mathrm{K}$. The results of the present model on a uniform mesh
with cell size $\Delta x=$ $0.01\text{\,}\mathrm{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%.

Figure 2: Code validation and solver verification. (a) One-dimensional
benchmark problem; time evolution of the mushy-zone thickness (blue: solidus;
red: liquidus) predicted by the present model implemented in Fluent and
OpenFOAM compared with available semi-analytical [51] and numerical [41]
results. (b) Two-dimensional benchmark problem; time evolution of solid–liquid
interface during isothermal phase-change with convection predicted by the
present model implemented in Fluent and OpenFOAM compared with numerical
results of Hannoun et al. [11].

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, $\Delta x_{\mathrm{i}}/\mathrm{W}=$
$4\text{\times}{10}^{-4}\text{\,}$ 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).

Figure 3: Comparisons of numerically determined melting-front positions in a
rectangular enclosure subject to natural convection using the present model
implemented in Fluent with experimental observations of isothermal melting of
lead using thermal neutron radiography [52].

###  5.2 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_{\mathrm{m}}$ = $302.78\text{\,}\mathrm{K}$, in a side-heated rectangular
enclosure with $T_{\mathrm{i}}=T_{\mathrm{c}}$ = $301.30\text{\,}\mathrm{K}$,
and $T_{\mathrm{h}}$ = $311.0\text{\,}\mathrm{K}$, which was experimentally
investigated by Gau and Viskanta [53]. This benchmark case has often been
considered for validation of phase-change 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 $\times$ 32, 105 $\times$
80, 210 $\times$ 160 and 420 $\times$ 320 computational grid cells and
$C=10^{4}$ and
$10^{8}\,$\text{\,}\mathrm{kg}\text{\,}{\mathrm{s}}^{-1}\text{\,}{\mathrm{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 solid–liquid 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
$\Delta\phi$ and $\Gamma$ (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\text{\times}{10}^{8}\text{\,}\mathrm{kg}\text{\,}{\mathrm{s}}^{-1}\text{\,}{\mathrm{m}}^{-3}$$
and
$C_{2}=$1\text{\times}{10}^{4}\text{\,}\mathrm{kg}\text{\,}{\mathrm{s}}^{-1}\text{\,}{\mathrm{m}}^{-3}$$
at $t=$150\text{\,}\mathrm{s}$$. Reducing the cell size decreases the
influence of the permeability coefficient on the results, with $\Gamma$
scaling approximately linearly with cell size, and $\Delta\phi$ 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_{\mathrm{L}}=0.0$) should be
large compared to the viscous term ($\nabla^{2}\vec{V}^{*}$) to suppress fluid
velocities (i.e., $C\,\mathrm{D}^{2}/\left(\mu\epsilon\right)\gg 1$). The
influence of the sink-term magnitude on predicted liquid fraction $\phi$ and
the ratio of the maximum velocity magnitude in solid regions to that in liquid
regions (i.e., $|\vec{V}_{\mathrm{solid}}|_{\mathrm{Max}}\ /\
|\vec{V}_{\mathrm{liquid}}|_{\mathrm{Max}}$) is shown in figure 6. Here,
$\mathrm{D}$ is chosen to be the size of the heated wall. For magnitudes of
$C\,\mathrm{D}^{2}/\left(\mu\epsilon\right)$ roughly larger than
$2.5\text{\times}{10}^{7}\text{\,}$, 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.

Figure 4: The influence of the computational cell size and the permeability
coefficient on predicted melting front position for the gallium melting
problem. Grid cell size: (a) $\Delta
x/\mathrm{D}\approx$3.33\text{\times}{10}^{-2}\text{\,}$$; (b) $\Delta
x/\mathrm{D}\approx$1.33\text{\times}{10}^{-2}\text{\,}$$; (c) $\Delta
x/\mathrm{D}\approx$6.67\text{\times}{10}^{-3}\text{\,}$$; and (d) $\Delta
x/\mathrm{D}\approx$3.33\text{\times}{10}^{-3}\text{\,}$$.

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].

Figure 5: The relative difference between the solid–liquid interface
morphologies ($\Gamma$ (defined in equation 13), circles in blue), and the
relative difference between the liquid fractions ($\Delta\phi$ (defined below
equation 12), squares in red) at $t=$150\text{\,}\mathrm{s}$$ when using
permeability coefficients
$C_{1}=$1\text{\times}{10}^{8}\text{\,}\mathrm{kg}\text{\,}{\mathrm{s}}^{-1}\text{\,}{\mathrm{m}}^{-3}$$
and
$C_{2}=$1\text{\times}{10}^{4}\text{\,}\mathrm{kg}\text{\,}{\mathrm{s}}^{-1}\text{\,}{\mathrm{m}}^{-3}$$
as a function of grid size ($\Delta x$). Larger values of $\Delta\phi$ and
$\Gamma$ indicate more sensitivity to the value of the permeability
coefficient $C$. Symbols: results obtained from numerical simulations; dashed
lines: curve-fit. Figure 6: The ratio of the maximum velocity magnitude in
solid regions to that in liquid regions
($|\vec{V}_{\mathrm{solid}}|_{\mathrm{Max}}\ /\
|\vec{V}_{\mathrm{liquid}}|_{\mathrm{Max}}$, blue circles) and the predicted
liquid fractions $\phi$ (red squares) for various values of the sink term
$\left(C\,\mathrm{D}^{2}/\left(\mu\epsilon\right)\right)$.

###  5.3 Non-Isothermal Phase-Change

####  5.3.1 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 $\Delta
T_{\mathrm{m}}=$50\text{\,}\mathrm{K}$$, figure 7 shows the sensitivity
parameter $\Delta\phi$, for short time (Fo = 0.12) and steady-state (Fo = 9.0)
obtained with
$C_{1}=$1\text{\times}{10}^{8}\text{\,}\mathrm{kg}\text{\,}{\mathrm{s}}^{-1}\text{\,}{\mathrm{m}}^{-3}$$
and
$C_{2}=$1\text{\times}{10}^{4}\text{\,}\mathrm{kg}\text{\,}{\mathrm{s}}^{-1}\text{\,}{\mathrm{m}}^{-3}$$.
For comparison, the same problem was also solved assuming isothermal phase
change ($\Delta T_{\mathrm{m}}=$0\text{\,}\mathrm{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 $\Delta
x_{\mathrm{i}}/\mathrm{D}=$4\text{\times}{10}^{-3}\text{\,}$$. Consequently,
this base grid size is used in the next sections. In addition, a four-level
dynamic solution-adaptive mesh refinement (i.e., $\Delta x_{\mathrm{N}}=\Delta
x_{\mathrm{i}}/2^{\mathrm{N}},\mathrm{N}=1,\cdots,4$) is applied to further
enhance the accuracy with which we capture the solid–liquid interface.

Figure 7: The relative difference between liquid fractions predicted with
different permeability coefficients (i.e.,
$C_{1}=$1\text{\times}{10}^{8}\text{\,}\mathrm{kg}\text{\,}{\mathrm{s}}^{-1}\text{\,}{\mathrm{m}}^{-3}$$
and
$C_{2}=$1\text{\times}{10}^{4}\text{\,}\mathrm{kg}\text{\,}{\mathrm{s}}^{-1}\text{\,}{\mathrm{m}}^{-3}$$)
on different grid densities. Diamonds (in blue): isothermal phase-change;
squares (in red): non-isothermal phase-change; unfilled symbols: Fo = 9.0
(steady-state condition); filled symbols: Fo = 0.12.

####  5.3.2 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_{\mathrm{L}}=1$) is shown in figure 8 for a fixed
wall temperature difference ($\Delta T_{\mathrm{w}}$), a fixed melting
temperature range ($\Delta T_{\mathrm{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 8: The influence of the permeability coefficient on predicted liquidus-
line ($f_{\mathrm{L}}$ = 1) position for different Rayleigh numbers ($\Delta
T_{\mathrm{m}}=$50\text{\,}\mathrm{K}$$; $\Delta
T_{\mathrm{w}}=$750\text{\,}\mathrm{K}$$; dashed lines:
$C=$1\text{\times}{10}^{8}\text{\,}\mathrm{kg}\text{\,}{\mathrm{s}}^{-1}\text{\,}{\mathrm{m}}^{-3}$$;
dotted lines:
$C=$1\text{\times}{10}^{4}\text{\,}\mathrm{kg}\text{\,}{\mathrm{s}}^{-1}\text{\,}{\mathrm{m}}^{-3}$$).
The results are given for Fo = 9, in which steady-state solutions are
achieved.

Figure 9 shows the effect of the permeability coefficient for a fixed melting
temperature range ($\Delta T_{\mathrm{m}}$) with different wall-temperature
differences ($\Delta T_{\mathrm{w}}$). Here, a higher value of $\Delta
T_{\mathrm{w}}$ (i.e., smaller $\theta$) leads to a smaller mushy zone
thickness, and consequently a reduced sensitivity to the permeability
coefficient. Similarly, for a fixed ($\Delta T_{\mathrm{w}}$), reducing
($\Delta T_{\mathrm{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 $\theta$), 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 ($\Delta T_{\mathrm{m}}$) and wall-temperature difference
($\Delta T_{\mathrm{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.

Figure 9: The influence of the permeability coefficient on predicted liquidus-
line ($f_{\mathrm{L}}$ = 1) position for different wall-temperature
differences $\Delta T_{\mathrm{w}}$, at $\Delta
T_{\mathrm{m}}=$50\text{\,}\mathrm{K}$$, and Ra = $10^{6}$. Dashed lines:
$C=$1\text{\times}{10}^{8}\text{\,}\mathrm{kg}\text{\,}{\mathrm{s}}^{-1}\text{\,}{\mathrm{m}}^{-3}$$;
dotted lines:
$C=$1\text{\times}{10}^{4}\text{\,}\mathrm{kg}\text{\,}{\mathrm{s}}^{-1}\text{\,}{\mathrm{m}}^{-3}$$.
The results are given for Fo = 9, in which steady-state solutions are
achieved. Figure 10: The influence of the permeability coefficient on
predicted liquidus-line ($f_{\mathrm{L}}$ = 1) position for different melting-
temperature ranges $\Delta T_{\mathrm{m}}$, at $\Delta
T_{\mathrm{w}}=$750\text{\,}\mathrm{K}$$, and Ra = $10^{6}$. Dashed lines:
$C=$1\text{\times}{10}^{8}\text{\,}\mathrm{kg}\text{\,}{\mathrm{s}}^{-1}\text{\,}{\mathrm{m}}^{-3}$$;
dotted lines:
$C=$1\text{\times}{10}^{4}\text{\,}\mathrm{kg}\text{\,}{\mathrm{s}}^{-1}\text{\,}{\mathrm{m}}^{-3}$$.
The results are given for Fo = 9, in which steady-state solutions are
achieved. Figure 11: The influence of the permeability coefficient on
predicted liquidus-line ($f_{\mathrm{L}}$ = 1) position for different values
of thermal diffusivity $\alpha$, at $\Delta
T_{\mathrm{m}}=$50\text{\,}\mathrm{K}$$, $\Delta
T_{\mathrm{w}}=$750\text{\,}\mathrm{K}$$ and Ra = $10^{6}$. Dashed lines:
$C=$1\text{\times}{10}^{8}\text{\,}\mathrm{kg}\text{\,}{\mathrm{s}}^{-1}\text{\,}{\mathrm{m}}^{-3}$$;
dotted lines:
$C=$1\text{\times}{10}^{4}\text{\,}\mathrm{kg}\text{\,}{\mathrm{s}}^{-1}\text{\,}{\mathrm{m}}^{-3}$$.
The results are given for Fo = 9, in which steady-state solutions are
achieved.

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
non-isothermal phase change with a thick mushy zone ($\theta$ = 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 ($\theta=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.

Figure 12: Time evolution of melting front positions ($f_{\mathrm{L}}$ = 1)
for isothermal phase-change (blue lines, $\theta$ = 0) and non-isothermal
phase-change (red lines: (a) $\Delta T_{\mathrm{m}}=$50\text{\,}\mathrm{K}$$,
$\Delta T_{\mathrm{w}}=$750\text{\,}\mathrm{K}$$, $\theta$ = 1/15; and (b)
$\Delta T_{\mathrm{m}}=$25\text{\,}\mathrm{K}$$, $\Delta
T_{\mathrm{w}}=$750\text{\,}\mathrm{K}$$, $\theta$ = 1/30) problems. Dashed
lines:
$C=$1\text{\times}{10}^{4}\text{\,}\mathrm{kg}\text{\,}{\mathrm{s}}^{-1}\text{\,}{\mathrm{m}}^{-3}$$;
dotted lines:
$C=$1\text{\times}{10}^{8}\text{\,}\mathrm{kg}\text{\,}{\mathrm{s}}^{-1}\text{\,}{\mathrm{m}}^{-3}$$
(Ra = $10^{6}$).

##  6 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_{\mathrm{t}}$) can be estimated as

$m_{\mathrm{t}}=\frac{\Delta T_{\mathrm{m}}}{\left|\nabla
T\right|}=\frac{\Delta
T_{\mathrm{m}}\cdot\mathscr{L}}{T_{\mathrm{h}}-T_{\mathrm{c}}},$ (14)

where $\mathscr{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, $\mathrm{Pe^{*}}$ is defined as follows:

$\mathrm{Pe^{*}}=\frac{|\vec{V}|\cdot
m_{\mathrm{t}}}{\alpha}=\frac{|\vec{V}|\cdot\Delta
T_{\mathrm{m}}\cdot\mathscr{L}}{\alpha\left(T_{\mathrm{h}}-T_{\mathrm{c}}\right)}.$
(15)

Non-zero values of $\mathrm{Pe^{*}}$ adjacent to the solid–liquid interface
indicate increased sensitivity to the permeability coefficient. For isothermal
phase-change, $\mathrm{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 $\Delta T_{\mathrm{m}}$ and a fixed
$\Delta T_{\mathrm{w}}=$750\text{\,}\mathrm{K}$$, and the results are shown in
figure 13. Higher values of $\mathrm{Pe^{*}}$ ($\approx\mathcal{O}(10)$) along
the melting front for the case with $\Delta
T_{\mathrm{m}}=$55\text{\,}\mathrm{K}$$ (figure 13c) compared to the case with
$\Delta T_{\mathrm{m}}=$10\text{\,}\mathrm{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 $\mathrm{Pe^{*}}$ along the interface are small (i.e.,
$\mathrm{Pe^{*}}\ll 1$), the results appear to be insensitive to the
permeability coefficient, while, for $\mathrm{Pe^{*}}\gg 1$, the results are
sensitive to $C$.

Figure 13: Melting front positions ($f_{\mathrm{L}}$ = 1) for the problem
described in section 2 ($\Delta T_{\mathrm{w}}=$750\text{\,}\mathrm{K}$$ and
Ra = $10^{6}$) at the steady-state condition (Fo = 9) and values of
$\mathrm{Pe^{*}}$ along the melting fronts for: (a) $\Delta
T_{\mathrm{m}}=$1\text{\,}\mathrm{K}$$; (b) $\Delta
T_{\mathrm{m}}=$10\text{\,}\mathrm{K}$$; and (c) $\Delta
T_{\mathrm{m}}=$55\text{\,}\mathrm{K}$$ predicted using different permeability
coefficients of
$C=10^{4}\,$\text{\,}\mathrm{kg}\text{\,}{\mathrm{s}}^{-1}\text{\,}{\mathrm{m}}^{-3}$$
(coloured lines), and
$10^{8}\,$\text{\,}\mathrm{kg}\text{\,}{\mathrm{s}}^{-1}\text{\,}{\mathrm{m}}^{-3}$$
(black lines). Values of $\mathrm{Pe^{*}}$ along the melting front are
$\mathcal{O}(0.1)$, $\mathcal{O}(1)$ and $\mathcal{O}(10)$ for $\Delta
T_{\mathrm{m}}=$1\text{\,}\mathrm{K}$$, $10\text{\,}\mathrm{K}$ and
$55\text{\,}\mathrm{K}$, respectively. The results are more sensitive to the
permeability coefficient for large $\mathrm{Pe^{*}}$. For the case with
$\Delta T_{\mathrm{m}}=$1\text{\,}\mathrm{K}$$ (a), predictions with different
permeability coefficient are identical and therefore the coloured line covers
the black line.

The general applicability of the proposed $\mathrm{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\text{\,}\mathrm{p}\mathrm{p}\mathrm{m}$ of sulphur was heated using a
stationary laser beam with a power of $5200\text{\,}\mathrm{W}$ and a top-hat
radius of $1.4\text{\,}\mathrm{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}\,$\text{\,}\mathrm{kg}\text{\,}{\mathrm{s}}^{-1}\text{\,}{\mathrm{m}}^{-3}$$.
To study the influence of the mushy-zone thickness on the sensitivity to $C$
through $\mathrm{Pe^{*}}$, the melting-temperature range ($\Delta
T_{\mathrm{m}}$) of the material is changed artificially to
$40\text{\,}\mathrm{K}$ and $200\text{\,}\mathrm{K}$.

Figure 14 shows the position of melting front ($f_{\mathrm{L}}$ = 1) at
$t=$5\text{\,}\mathrm{s}$$, as well as the value of $\mathrm{Pe^{*}}$ along it
when using different permeability coefficients and melting temperature ranges.
The laser-beam diameter is chosen here as the characteristic length scale
$\mathscr{L}$ to calculate $\mathrm{Pe^{*}}$. Higher values of
$\mathrm{Pe^{*}}$ are found for larger $\Delta T_{\mathrm{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.

Figure 14: Laser spot melting of Fe-S binary alloys with
$20\text{\,}\mathrm{p}\mathrm{p}\mathrm{m}$ sulphur content. Melting front
positions ($f_{\mathrm{L}}$ = 1) at $t=$5\text{\,}\mathrm{s}$$ and values of
$\mathrm{Pe^{*}}$ along the melting fronts, predicted using different
permeability coefficients of
$C=10^{6}\,$\text{\,}\mathrm{kg}\text{\,}{\mathrm{s}}^{-1}\text{\,}{\mathrm{m}}^{-3}$$
(coloured lines), and
$10^{10}\,$\text{\,}\mathrm{kg}\text{\,}{\mathrm{s}}^{-1}\text{\,}{\mathrm{m}}^{-3}$$
(black lines): (a) $\Delta T_{\mathrm{m}}=$200\text{\,}\mathrm{K}$$; and (b)
$\Delta T_{\mathrm{m}}=$40\text{\,}\mathrm{K}$$. Orange symbols show the
experimental observation of the melting front position reported in [60].

##  7 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 enthalpy-porosity 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 $\mathrm{Pe^{*}}\gg 1$ adjacent to the solid–liquid
interface indicate a strong sensitivity and $\mathrm{Pe^{*}}\ll 1$ indicates
insensitivity to the permeability coefficient $C$.

## Author Contributions

Conceptualisation, A.E., C.R.K. and I.M.R.; methodology, A.E.; software, A.E.;
validation, A.E.; formal analysis, A.E.; investigation, A.E.; resources, A.E.,
C.R.K, and I.M.R.; data curation, A.E.; writing—original draft preparation,
A.E.; writing—review and editing, A.E., C.R.K., and I.M.R; visualisation,
A.E.; supervision, C.R.K. and I.M.R.; project administration, A.E. and I.M.R.;
and funding acquisition, I.M.R.

## Acknowledgement

This research was carried out under project number F31.7.13504 in the
framework of the Partnership Program of the Materials innovation institute M2i
(www.m2i.nl) and the Foundation for Fundamental Research on Matter (FOM)
(www.fom.nl), which is part of the Netherlands Organisation for Scientific
Research (www.nwo.nl). The authors would like to thank the industrial partner
in this project “Allseas Engineering B.V.” for the financial support.

## References

  * [1] J. Crank, Free and moving boundary problems, Clarendon Press, Oxford, UK, 1984.
  * [2] W. Shyy, Multi-scale computational heat transfer with moving solidification boundaries, International Journal of Heat and Fluid Flow 23 (3) (2002) 278–287.  doi:10.1016/s0142-727x(02)00175-3.
  * [3] M. Rappaz, Modelling of microstructure formation in solidification processes, International Materials Reviews 34 (1) (1989) 93–124.  doi:10.1179/imr.1989.34.1.93.
  * [4] Y. Dutil, D. R. Rousse, N. B. Salah, S. Lassue, L. Zalewski, A review on phase-change materials: Mathematical modeling and simulations, Renewable and Sustainable Energy Reviews 15 (1) (2011) 112–130.  doi:10.1016/j.rser.2010.06.011.
  * [5] S. Verma, A. Dewan, Solidification modeling: Evolution, benchmarks, trends in handling turbulence, and future directions, Metallurgical and Materials Transactions B 45 (4) (2014) 1456–1471.  doi:10.1007/s11663-014-0039-6.
  * [6] M. A. Jaafar, D. R. Rousse, S. Gibout, J.-P. Bédécarrats, A review of dendritic growth during solidification: Mathematical modeling and numerical simulations, Renewable and Sustainable Energy Reviews 74 (2017) 1064–1079.  doi:10.1016/j.rser.2017.02.050.
  * [7] B. Basu, A. W. Date, Numerical modelling of melting and solidification problems—a review, Sadhana 13 (3) (1988) 169–213.  doi:10.1007/bf02812200.
  * [8] M. Lacroix, V. R. Voller, Finite difference solutions of solidification phase change problems: transformed versus fixed grids, Numerical Heat Transfer, Part B: Fundamentals 17 (1) (1990) 25–41.  doi:10.1080/10407799008961731.
  * [9] V. R. Voller, Numerical Methods for Phase-Change Problems, Wiley-Blackwell, 2009, Ch. 19, pp. 593–622.  doi:10.1002/9780470172599.ch19.
  * [10] J. Mencinger, Numerical simulation of melting in two-dimensional cavity using adaptive grid, Journal of Computational Physics 198 (1) (2004) 243–264.  doi:10.1016/j.jcp.2004.01.006.
  * [11] N. Hannoun, V. Alexiades, T. Z. Mai, A reference solution for phase change with convection, International Journal for Numerical Methods in Fluids 48 (11) (2005) 1283–1308.  doi:10.1002/fld.979.
  * [12] C. Lan, C. Liu, C. Hsu, An adaptive finite volume method for incompressible heat flow problems in solidification, Journal of Computational Physics 178 (2) (2002) 464–497.  doi:10.1006/jcph.2002.7037.
  * [13] M. G. Worster, Convection in mushy layers, Annual Review of Fluid Mechanics 29 (1) (1997) 91–122.  doi:10.1146/annurev.fluid.29.1.91.
  * [14] P. L. Quéré, D. Gobin, A note on possible flow instabilities in melting from the side, International Journal of Thermal Sciences 38 (7) (1999) 595–600.  doi:10.1016/s0035-3159(99)80039-7.
  * [15] K. Morgan, A numerical analysis of freezing and melting with convection, Computer Methods in Applied Mechanics and Engineering 28 (3) (1981) 275–284.  doi:10.1016/0045-7825(81)90002-5.
  * [16] V. Voller, C. Prakash, A fixed grid numerical modelling methodology for convection-diffusion mushy region phase-change problems, International Journal of Heat and Mass Transfer 30 (8) (1987) 1709–1719.  doi:10.1016/0017-9310(87)90317-6.
  * [17] V. R. Voller, M. Cross, N. C. Markatos, An enthalpy method for convection/diffusion phase change, International Journal for Numerical Methods in Engineering 24 (1) (1987) 271–284.  doi:10.1002/nme.1620240119.
  * [18] H. C. Brinkman, A calculation of the viscous force exerted by a flowing fluid on a dense swarm of particles, Flow, Turbulence and Combustion 1 (1).  doi:10.1007/bf02120313.
  * [19] D. R. Poirier, Permeability for flow of interdendritic liquid in columnar-dendritic alloys, Metallurgical Transactions B 18 (1) (1987) 245–255.  doi:10.1007/bf02658450.
  * [20] A. K. Singh, R. Pardeshi, B. Basu, Modelling of convection during solidification of metal and alloys, Sadhana 26 (1-2) (2001) 139–162.  doi:10.1007/bf02728483.
  * [21] M. Fadl, P. C. Eames, Numerical investigation of the influence of mushy zone parameter amush on heat transfer characteristics in vertically and horizontally oriented thermal energy storage systems, Applied Thermal Engineering 151 (2019) 90–99.  doi:10.1016/j.applthermaleng.2019.01.102.
  * [22] Y. Hong, W.-B. Ye, J. Du, S.-M. Huang, Solid-liquid phase-change thermal storage and release behaviors in a rectangular cavity under the impacts of mushy region and low gravity, International Journal of Heat and Mass Transfer 130 (2019) 1120–1132.  doi:10.1016/j.ijheatmasstransfer.2018.11.024.
  * [23] R. Rai, G. G. Roy, T. DebRoy, A computationally efficient model of convective heat transfer and solidification characteristics during keyhole mode laser welding, Journal of Applied Physics 101 (5) (2007) 054909.  doi:10.1063/1.2537587.
  * [24] Y. Zheng, Q. Li, Z. Zheng, J. Zhu, P. Cao, Modeling the impact, flattening and solidification of a molten droplet on a solid substrate during plasma spraying, Applied Surface Science 317 (2014) 526–533.  doi:10.1016/j.apsusc.2014.08.032.
  * [25] X.-H. Yang, S.-C. Tan, J. Liu, Numerical investigation of the phase change process of low melting point metal, International Journal of Heat and Mass Transfer 100 (2016) 899–907.  doi:10.1016/j.ijheatmasstransfer.2016.04.109.
  * [26] T. Kousksou, M. Mahdaoui, A. Ahmed, A. A. Msaad, Melting over a wavy surface in a rectangular cavity heated from below, Energy 64 (2014) 212–219.  doi:10.1016/j.energy.2013.11.033.
  * [27] R. Karami, B. Kamkari, Investigation of the effect of inclination angle on the melting enhancement of phase change material in finned latent heat thermal storage units, Applied Thermal Engineering 146 (2019) 45–60.  doi:10.1016/j.applthermaleng.2018.09.105.
  * [28] C. Pan, J. Charles, N. Vermaak, C. Romero, S. Neti, Y. Zheng, C.-H. Chen, R. Bonner, Experimental, numerical and analytic study of unconstrained melting in a vertical cylinder with a focus on mushy region effects, International Journal of Heat and Mass Transfer 124 (2018) 1015–1024.  doi:10.1016/j.ijheatmasstransfer.2018.04.009.
  * [29] S. Arena, E. Casti, J. Gasia, L. F. Cabeza, G. Cau, Numerical simulation of a finned-tube LHTES system: influence of the mushy zone constant on the phase change behaviour, Energy Procedia 126 (2017) 517–524.  doi:10.1016/j.egypro.2017.08.237.
  * [30] M. Prieto, B. González, Fluid flow and heat transfer in PCM panels arranged vertically and horizontally for application in heating systems, Renewable Energy 97 (2016) 331–343.  doi:10.1016/j.renene.2016.05.089.
  * [31] A. C. Kheirabadi, D. Groulx, Simulating phase change heat transfer using COMSOL and Fluent: Effect of the mushy-zone constant, Computational Thermal Sciences: An International Journal 7 (5-6) (2015) 427–440.  doi:10.1615/computthermalscien.2016014279.
  * [32] S. Hosseinizadeh, A. R. Darzi, F. Tan, J. Khodadadi, Unconstrained melting inside a sphere, International Journal of Thermal Sciences 63 (2013) 55–64.  doi:10.1016/j.ijthermalsci.2012.07.012.
  * [33] H. Shmueli, G. Ziskind, R. Letan, Melting in a vertical cylindrical tube: Numerical investigation and comparison with experiments, International Journal of Heat and Mass Transfer 53 (19-20) (2010) 4082–4091.  doi:10.1016/j.ijheatmasstransfer.2010.05.028.
  * [34] M. Kumar, D. J. Krishna, Influence of mushy zone constant on thermohydraulics of a PCM, Energy Procedia 109 (2017) 314–321.  doi:10.1016/j.egypro.2017.03.074.
  * [35] H. Sattari, A. Mohebbi, M. Afsahi, A. A. Yancheshme, CFD simulation of melting process of phase change materials (PCMs) in a spherical capsule, International Journal of Refrigeration 73 (2017) 209–218.  doi:10.1016/j.ijrefrig.2016.09.007.
  * [36] Z. Hu, A. Li, R. Gao, H. Yin, Effect of the length ratio on thermal energy storage in wedge-shaped enclosures, Journal of Thermal Analysis and Calorimetry 117 (2) (2014) 807–816.  doi:10.1007/s10973-014-3843-y.
  * [37] J. Vogel, J. Felbinger, M. Johnson, Natural convection in high temperature flat plate latent heat thermal energy storage systems, Applied Energy 184 (2016) 184–196.  doi:10.1016/j.apenergy.2016.10.001.
  * [38] M. Hameter, H. Walter, Influence of the mushy zone constant on the numerical simulation of the melting and solidification process of phase change materials, in: Computer Aided Chemical Engineering, Elsevier, 2016, pp. 439–444.  doi:10.1016/b978-0-444-63428-3.50078-3.
  * [39] D. J. Tritton, Physical Fluid Dynamics, 1st Edition, Springer Netherlands, 1977\.  doi:10.1007/978-94-009-9992-3.
  * [40] V. R. Voller, C. R. Swaminathan, B. G. Thomas, Fixed grid techniques for phase change problems: A review, International Journal for Numerical Methods in Engineering 30 (4) (1990) 875–898.  doi:10.1002/nme.1620300419.
  * [41] V. R. Voller, C. R. Swaminathan, General source-based method for solidification phase change, Numerical Heat Transfer, Part B: Fundamentals 19 (2) (1991) 175–189.  doi:10.1080/10407799108944962.
  * [42] C. R. Swaminathan, V. R. Voller, A general enthalpy method for modeling solidification processes, Metallurgical Transactions B 23 (5) (1992) 651–664.  doi:10.1007/bf02649725.
  * [43] I. L. Ferreira, V. R. Voller, B. Nestler, A. Garcia, Two-dimensional numerical model for the analysis of macrosegregation during solidification, Computational Materials Science 46 (2) (2009) 358–366.  doi:10.1016/j.commatsci.2009.03.020.
  * [44] M. Faraji, H. E. Qarnia, Numerical study of melting in an enclosure with discrete protruding heat sources, Applied Mathematical Modelling 34 (5) (2010) 1258–1275.  doi:10.1016/j.apm.2009.08.012.
  * [45] S. Bouabdallah, R. Bessaih, Effect of magnetic field on 3D flow and heat transfer during solidification from a melt, International Journal of Heat and Fluid Flow 37 (2012) 154–166.  doi:10.1016/j.ijheatfluidflow.2012.07.002.
  * [46] M. Mahdaoui, T. Kousksou, S. Blancher, A. A. Msaad, T. E. Rhafiki, M. Mouqallid, A numerical analysis of solid–liquid phase change heat transfer around a horizontal cylinder, Applied Mathematical Modelling 38 (3) (2014) 1101–1110.  doi:10.1016/j.apm.2013.08.002.
  * [47] R. Y. Farsani, A. Raisi, A. A. Nadooshan, S. Vanapalli, Does nanoparticles dispersed in a phase change material improve melting characteristics?, International Communications in Heat and Mass Transfer 89 (2017) 219–229.  doi:10.1016/j.icheatmasstransfer.2017.10.006.
  * [48] J. Dannenhoffer, III, J. Baron, Grid adaptation for the 2-D Euler equations, in: 23rd Aerospace Sciences Meeting, American Institute of Aeronautics and Astronautics, 1985.  doi:10.2514/6.1985-484.
  * [49] R. Issa, Solution of the implicitly discretised fluid flow equations by operator-splitting, Journal of Computational Physics 62 (1) (1986) 40–65.  doi:10.1016/0021-9991(86)90099-9.
  * [50] S. V. Patankar, Numerical Heat Transfer and Fluid Flow, 1st Edition, Taylor & Francis Inc, 1980.
  * [51] V. Voller, Development and application of a heat balance integral method for analysis of metallurgical solidification, Applied Mathematical Modelling 13 (1) (1989) 3–11.  doi:10.1016/0307-904x(89)90191-1.
  * [52] L. Kumar, B. Manjunath, R. Patel, S. Markandeya, R. Agrawal, A. Agrawal, Y. Kashyap, P. Sarkar, A. Sinha, K. Iyer, S. Prabhu, Experimental investigations on melting of lead in a cuboid with constant heat flux boundary condition using thermal neutron radiography, International Journal of Thermal Sciences 61 (2012) 15–27.  doi:10.1016/j.ijthermalsci.2012.06.014.
  * [53] C. Gau, R. Viskanta, Melting and solidification of a pure metal on a vertical wall, Journal of Heat Transfer 108 (1) (1986) 174.  doi:10.1115/1.3246884.
  * [54] N. Hannoun, V. Alexiades, T. Z. Mai, Resolving the controversy over tin and gallium melting in a rectangular cavity heated from the side, Numerical Heat Transfer, Part B: Fundamentals 44 (3) (2003) 253–276.  doi:10.1080/713836378.
  * [55] Y. Lee, S. A. Korpela, Multicellular natural convection in a vertical slot, Journal of Fluid Mechanics 126 (-1) (1983) 91.  doi:10.1017/s0022112083000063.
  * [56] J. A. Dantzig, Modelling liquid-solid phase changes with melt convection, International Journal for Numerical Methods in Engineering 28 (8) (1989) 1769–1785.  doi:10.1002/nme.1620280805.
  * [57] M. Cerimele, D. Mansutti, F. Pistella, Numerical modelling of liquid/solid phase transitions: Analysis of a gallium melting test, Computers & Fluids 31 (4-7) (2002) 437–451.  doi:10.1016/s0045-7930(01)00062-7.
  * [58] P. W. Schroeder, G. Lube, Stabilised dG-FEM for incompressible natural convection flows with boundary and moving interior layers on non-adapted meshes, Journal of Computational Physics 335 (2017) 760–779.  doi:10.1016/j.jcp.2017.01.055.
  * [59] J. Vogel, A. Thess, Validation of a numerical model with a benchmark experiment for melting governed by natural convection in latent thermal energy storage, Applied Thermal Engineering 148 (2019) 147–159.  doi:10.1016/j.applthermaleng.2018.11.032.
  * [60] W. Pitscheneder, T. Debroy, K. Mundra, R. Ebner, Role of sulfur and processing variables on the temporal evolution of weld pool geometry during multikilowatt laser beam welding of steels, Welding Journal 75 (3) (1996) 71–80.
  * [61] J. Yan, W. Yan, S. Lin, G. Wagner, A fully coupled finite element formulation for liquid–solid–gas thermo-fluid flow with melting and solidification, Computer Methods in Applied Mechanics and Engineering 336 (2018) 444–470.  doi:10.1016/j.cma.2018.03.017.
  * [62] Z. Saldi, A. Kidess, S. Kenjereš, C. Zhao, I. Richardson, C. Kleijn, Effect of enhanced heat and mass transport and flow reversal during cool down on weld pool shapes in laser spot welding of steel, International Journal of Heat and Mass Transfer 66 (2013) 879–888.  doi:10.1016/j.ijheatmasstransfer.2013.07.085.
  * [63] P. Sahoo, T. Debroy, M. J. McNallan, Surface tension of binary metal—surface active solute systems under conditions relevant to welding metallurgy, Metallurgical Transactions B 19 (3) (1988) 483–491.  doi:10.1007/bf02657748.


11institutetext: College of Computer Science, Chongqing University,
Chongqing 400044, China
11email: $\\{$majunhua<EMAIL_ADDRESS><EMAIL_ADDRESS>School of Information Technology &
Electrical Engineering,
The University of Queensland, St Lucia Qld 4072, Australia
22email<EMAIL_ADDRESS>
Junhua Ma, Jiajun Li, Xueming Li and Xu Li

# PathSAGE: Spatial Graph Attention Neural Networks With Random Path Sampling

Junhua Ma 11 Jiajun Li 22 Xueming Li 11 Xu Li 11

###### Abstract

Graph Convolutional Networks (GCNs) achieve great success in non-Euclidean
structure data processing recently. In existing studies, deeper layers are
used in CCNs to extract deeper features of Euclidean structure data. However,
for non-Euclidean structure data, too deep GCNs will confront with problems
like ”neighbor explosion” and ”over-smoothing”, it also cannot be applied to
large datasets. To address these problems, we propose a model called PathSAGE,
which can learn high-order topological information and improve the model’s
performance by expanding the receptive field. The model randomly samples paths
starting from the central node and aggregates them by Transformer encoder.
PathSAGE has only one layer of structure to aggregate nodes which avoid those
problems above. The results of evaluation shows that our model achieves
comparable performance with the state-of-the-art models in inductive learning
tasks.

###### Keywords:

Neural network models, Path, GCNs, Transformer, Random

##  1 Introduction

Convolutional Neural Networks (CNNs) have been successfully used in various
tasks with Euclidean structure data in recent years. For non-Euclidean
structure datasets, graph convolutional networks (GCNs) use the same idea to
extract the topological structure information.

[1] firstly proposed two GCN models with spectral and spatial construction
respectively. The spectral model using a Laplacian matrix to aggregate
neighborhood information of each node in a graph. The spatial model partition
graph into clustering and update them by aggregating function. In order to
extract deeper features, model based on CNNs usually deepen the model’s
layers. While in GCNs, Deepening layers will cause a lot of problems. In
spectral construction method, too many layers lead to “over smoothing”[13,
25]: the features of nodes in graph will tend to be the same. In spatial
construction method, it will cause exponential growth of the number of sampled
neighbor nodes, called “neighbor explosion”. Node sampling and layer
sampling[7, 3, 6, 9, 20] were proposed to handle this problem, but due to
incomplete sampling, the inaccuracy of nodes’ representation accumulates
errors between layers.

To address these two problems, in this paper, we propose a model called
PathSAGE, which can learn high-order topological information by expanding the
receptive field. Firstly, We design a path sampling technique based on random
walk to sample paths starting from central nodes with different lengths, then
the sequences of paths are fed into Transformer encoder[17], which can extract
the semantic and distance information in sequence effectively. As shown in
Fig. 1, We view the sequences of paths from the tail nodes to the central node
as the central node’ neighbors. Secondly, following this idea, we take the
average of paths of the same length as the representation of the central node
in this level of reception field. Finally, after concatenating the aggregated
features of paths with different lengths, the final representation of the
central node are used for downstream tasks.

For the two problems mentioned earlier, on the one hand, the aggregation of
the central node only perform once in training of a sample, which never cased
“over-smoothing”. On the other hand, all the paths were sampled with a fixed
number for representing the central node in our model, instead of recursively
sampling exponentially growing neighbor nodes. And each path only contributes
to the central node, we do not need to calculate and store the temporary
representations of nodes from the middle layer. Furthermore, it prevents the
error propagation caused by incomplete sampling.

Our contribution can be summarized in three points:

  * •

We utilize the path sampling to take place of the node sampling to avoid error
accumulation caused by incomplete node sampling.

  * •

We propose and evaluate our model Path-SAGE to solve the existing “neighour
explosion” and “over-smoothing” problems. the model can capture richer and
more diverse patterns around the central node with only one layer of structure
to the central node.

  * •

We evaluate our model on three inductive learning tasks, and it reaches the
state-of-the-art performance on two of them. We analyze the attention weights
of Transformer encoder and detect some patterns in the attention mechanism,
which can further illustrate how the model works.

##  2 Related work

GNNs model was initially proposed by [1], and the convolution operation in the
traditional Euclidean structure was introduced into the graph network with the
non-Euclidean structure in this article. They [1] divided GNN into two
construction methods: spectral construction and spatial construction.
Subsequently, many studies are carried out around these two aspects.

In spectral construction, [5] used Chebyshev polynomials with learnable
parameters to approximate a smooth filter in the spectral domain, which
improves computation efficiency significantly. [11] further reduced the
computational cost through local first-order approximation. In spatial
construction, MoNet was proposed in [15], developing the GCN model by defining
a different weighted sum of the convolution operation, using the weighted sum
of the nodes as the central node feature instead of the average value. [7]
attempted various aggregator to gather the features of neighbor nodes. [18,
23] defined the convolution operation with a self-attention mechanism between
the central node and neighbor nodes, [14] brought the LSTM(Long Short Term
Memory networks)[8] from NLP (Natural Language Processing) to GNN, and built
an adaptive depth structure by applying the memory gates of LSTM. [19] sampled
the shortest path based on their attention score to the central node and
combined them in a mixed way. PPNP and APPNP[12] used the relationship between
GCN and PageRank to derive an improved communication scheme based on
personalized PageRank.

With the scale of graph data increasing, the full-batch algorithm is no longer
applicable, and the mini-batch algorithm using stochastic gradient descent is
applied. In recent years, some studies based on different sampling strategies
have been proposed. [7] tried random neighbor node sampling for the first time
to limit large amounts of nodes caused by recursive node sampling. Layer-wise
sampling techniques was applied in [3, 6, 9, 20], which only consider fixed
number of neighbors in the graph to avoid “neighbor explosion”. [3, 2] used a
control variate-based algorithm that can achieve good convergence by reducing
the approximate variance. Besides, subgraph sampling was first introduced by
[21], which used a probability distribution based on degree. [4] performed
clustering decomposition on the graph before the training phase and randomly
sampled a subgraph as a training batch at every step. [22] sampled subgraph
during training based on node and edge importance. All these sampling
techniques are based on GCNs, how to extract deeper topology information from
large graph datasets is still a problem.

##  3 Proposed method

In this section, we present the PathSAGE. Firstly, the sampling algorithm is
introduced in 3.1. Secondly, we detail the aggregator in 3.2. Finally, we
discuss the difference between our model and the related models.

Input: garph $G(V,E)$ central node $c$

sample depth $s$

sample num each length $L=\\{n_{1},n_{2},…,n_{s}\\}$}

Output: path sequcences with different length

$\\{\bm{P}_{1},\bm{P}_{2},...,\bm{P}_{s}\\}$

1 Ramdom Path Sampling(_$G$ , $c$, $s$, $L$_)

2foreach  _$l=0\to s$_ do

3 $i\leftarrow 0$;

4 while  _$i <n_{l}$_ do

5 $u\leftarrow c$;

6 $P\leftarrow\\{u\\}$;

7 for  _$j=0\to l$_ do

8 $u\leftarrow$ Node randomly selected from $u$’s neighbors;

9 $P\leftarrow P\cup\\{u\\}$

10 end for

11

12 $\bm{P}_{l}\leftarrow\bm{P}_{l}\cup\\{P\\}$;

13 $i\leftarrow i+1$;

14 end while

15

16 end foreach

return $\\{\bm{P}_{1},\bm{P}_{2},...,\bm{P}_{s}\\}$

Algorithm 1 Random Path Sampling

###  3.1 Path Sampling

Except deepening model, another way to expand receptive field of CNNs is to
increase the size of convolution kernel. Following this idea, we sample the
node sequences starting from central node and regard these paths as the
context of the central node. Therefore, the receptive field can be expanded by
extending these paths.

We use a straightforward random sampling algorithm based on random walk, shown
in Algorithm 1: for a central node, a sampling starts from it and randomly
selects a neighbor node each time until reaching the preset path length, and
multiple path sequences for the corresponding central node of various lengths
can be obtained in this way, constituting a training sample.

Figure 1: Mechanism of sampling and aggregation. (a) and (b) are two different
possible training samples with the same central node. Paths with same colors
are with same lengths same length and share the same aggregators.  Figure 2:
Structure of first aggregation in aggregator. A specific sampled path’s
features adds the position embedding, through Transformer encoder layers to
fuse the information.

###  3.2 Aggregator

There are two aggregations in the aggregator: the first one aggregate path
sequences; the second one aggregate different paths as the final
representation of the central node. For the first aggregation, We formulate
each path sequence as:

${\bm{P}}_{ij}\;=\;\\{{\vec{a}_{ij}^{1},\;\vec{a}_{ij}^{2},\;\vec{a}_{ij}^{3},\;...,\;\vec{a}_{ij}^{j}\;}\\}$
(1)

To utilize the position information of path sequence, we define positions of
nodes in paths as their distances to the central nodes. As same as
Transformer, we add nodes’ features and their positional vectors together as
the input of the structure:

$\bm{\tilde{P}}_{ij}\;=\;\\{{a_{ij}^{1}\;+\;pos\\_emb\left(1\right),...,\;a_{ij}^{j}\;+\;pos\\_emb\left(j\right)}\\}$
(2)

For the $pos\\_emb(\cdot)$, we use $sine$ and $cosine$ functions of different
frequencies to generate the positional vectors. In each dimension, the
positional embedding is:

$pos\\_emb(p)_{2i}=sin(\frac{{p}}{{{{10000}^{\frac{{2i}}{{{d}}}}}}})$ (3)

$pos\\_emb(p)_{2i+1}=cos(\frac{{p}}{{{{10000}^{\frac{{2i}}{{d}}}}}})$ (4)

Where $p$ is the position in the sequences, $2i$ and $2i+1$ mean the odd and
even dimensions of position embedding, and $d$ is the number of features.

After that, we apply the Transformer encoder on each path:

$\bm{\tilde{P}}_{ij}^{k}=transformer\\_bloc{k^{k}}(\bm{\tilde{P}}_{ij}^{k-1})$
(5)

Where $k$ means $k$-th Transformer encoder layer. Noted that, in a $m$-layer
Transformer encoder, the output of the last layer $\bm{\tilde{P}}_{ij}^{m}$ is
a sequence of features. We only take the output at position 0 as the final
representation of the path sequence.

$\tilde{P}_{ij}^{\prime}=\left[\bm{\tilde{P}}_{ij}^{m}\right]_{0}$ (6)

Following these equations (2)-(6), we can obtain representations of paths with
different length to the central node, as shown in Fig. 1. In the second
aggregation, we apply an average pooling layer to aggregate the central node
paths. Then we concatenate all the path representation and apply a feed-
forward layer with non-linearity to fuse these features. The final output of a
central node $C^{\prime}$ is computed as following:

$\begin{array}[]{c}C=\;concat\left({{C_{1}},\;...,\;{C_{s}}}\right)\\\ \\\
where\;{C_{i}}\;=\;Average\left(\tilde{P}_{i1}^{\prime},\;...,\;\tilde{P}_{in}^{\prime}\right)\end{array}$
(7)

${C^{\prime}}\;=\;max\left({0,\;C{\bm{W}_{1}}\;+\,{b_{1}}}\right){\bm{W}_{2}}\;+\;{b_{2}}$
(8)

where $s$ denotes the sample depth, $n$ denote the number of paths sampled in

###  3.3 Comparisons to related work

  * •

We introduce a sequence transduction model into our structure, but it is
distinct from related work based on these models. LSTM[8] was also used in
GraphSAGE[22] to aggregate node’s features, it is very sensitive to the order
of the input sequence. In contrast, neighbor nodes are disordered, and the
authors have rectified it by consistently feeding randomly-ordered sequences
to the LSTM. In our model, paths to the central node already have orders,
which are naturally suitable for sequence processing model.

  * •

The difference between our attention mechanism and GAT is that we collect all
the mutual attention information of each node in the path. The output of our
model is integrated information of entire sequence, instead of only attend to
the central node.

  * •

Path sampling method is also used in SPAGAN[19]. But each sampled path in
SPAGAN is the shortest one to the central node, this sample technique will
leads to a high computational overhead and limits the ability of the model to
be applied to large graph datasets. By contrast, we sample paths randomly,
which greatly save the computation overhead. Besides, in our sample algorithm,
the same node may have a different path to the central node, which may help
the model acquire more diverse patterns around the central node while saving
the computational overhead.

##  4 Experiments

In this section, we introduce the datasets and the experiment setting in 4.1
and 4.2 respectively. we present the results of evaluation in 4.3.

Table 1: Summary of inductive learning tasks’ statistics. | Reddit | Yelp | Flickr
---|---|---|---
type | single-label | multi-label | multi-label
#Node | 232,965 | 716,847 | 89,250
#Edges | 11,606,919 | 6,977,410 | 899,756
#Features | 602 | 300 | 500
#Classes | 41 | 100 | 7
Train / Val / Test |  66% / 10% / 24% |  75% / 10% / 15% |  50% / 25% / 25%

###  4.1 Dataset

Three large-scale datasets are used to evaluate our model: 1) Reddit[7]: a
collection of monthly user interaction networks from the year 2014 for 2046
sub-reddit communities from Reddit that connect users when one has replied to
the other, 2) Flickr: a social network built by forming links between images
sharing common metadata from Flickr. 3) Yelp[22]: a network that links up
businesses with the interaction of its customers. Reddit is a multiclass node-
wise classification task; Flickr and Yelp are multilabel classification tasks.
The detail statistics of these datasets are shown in Tab 1.

Table 2: Perforamnce on inductive learning tasks(Micro-F1). Model | Reddit | Yelp | Flickr
---|---|---|---
GCN | 0.933$\pm$0.000 | 0.378$\pm$0.001 | 0.492$\pm$0.003
GraphSAGE | 0.953$\pm$0.001 | 0.634$\pm$0.006 | 0.501$\pm$0.013
FastGCN | 0.924$\pm$0.001 | 0.265$\pm$0.053 | 0.504$\pm$0.001
S-GCN | 0.964$\pm$0.001 | 0.640$\pm$0.002 | 0.482$\pm$0.003
AS-GCN | 0.958$\pm$0.001 | - | 0.504$\pm$0.002
ClusterGCN | 0.954$\pm$0.001 | 0.609$\pm$0.005 | 0.481$\pm$0.005
GraphSAINT | 0.966$\pm$0.001 | 0.653$\pm$0.003 | 0.511$\pm$0.001
ours | 0.969$\pm$0.002 | 0.642$\pm$0.005 | 0.511$\pm$0.003

###  4.2 Experiment setup

We build our model on Pytorch framework[16] and construct the Transformer
encoder based on UER[24]. For all tasks, We train the model with Adam SGD
optimizer[10] with learning rate 1e-3 and a learning rate scheduler with 0.1
warmup ratio. We use two layers Transformer, each of which has 8 attention
heads to gather the features of paths. The batch size is 32. The dropout
layers are applied between each sub-layer in Transformer layer. We use
different dropout rates in the output layer and Transformer encoder, which are
0.3 and 0.1 respectively. We set the sampling length of the path ranging from
1 to 8 (depth $s=8$), and the number of paths sampled in each length are [5,
5, 5, 5, 5, 10, 10, 10]. For multi-label classification task Flickr and Yelp,
the final output is obtained through a $sigmoid$ activation, and in Reddit,
the final output is obtained through a $softmax$ activation. The hidden
dimension is 128 for Reddit and Flickr, 512 for Yelp.

###  4.3 Result

We compare our model with seven state-of-the-art model: GCN[11], GraphSAGE[7],
FastGCN[3], S-GCN[2], AS-GCN[9], ClusterGCN[4], GraphSAINT[22]. GraphSAGE uses
a random node sampling and LSTM aggregator. FastGCN and S-GCN use a control
variate-based algorithm that can achieve good convergence by reducing the
sampling variance. ClusterGCN performs clustering decomposition on the graph
before the training phase and randomly sampling a subgraph as a training-batch
at every step. GraphSAINT samples subgraph while training with several
sampling methods based on node and edge importance. The evaluation on these
tasks uses the $Micro-F1$ metric, and report the mean and confidence interval
of the metrics by five runs.

The results of inductive learning experiments are shown in Tab 2. As we can
see from the table, For Reddit, our model outperforms all the baseline models.
For Flickr, we achieve a comparable F1-score with the top-performing model.
For Yelp, we surpass most GCN models, second only to GraphSAINT. One
hypothesis to explain the difference of the results is: Reddit and Flickr have
more training samples and number of nodes features, which makes the attention
mechanism have enough data to capture the relationships between nodes.

Figure 3: Detection of attention weights in Transformer that node have the
same label in the path sequence may receive higher weights at some attention
heads. The heatmap in the upper-left corner is the weight of position 4 of
different layers attend to position 2 and attention heads, the right one is
the complete heatmap of the selected head. The lower-left picture shows the
labels of the nodes in the sequence.

##  5 Attention Analysis

We observed the attention weight of the trained model in the PubMed test set.
At some attention heads, we find that nodes with the same labels get very high
attention scores to each other. We visualize an example in Fig 3. Central node
0 and node 4 have the same label and receive extremely high attention scores
on second attention heads of the first layer. By observing the attention
scores for the whole sequence, we can see that this score also occupies a
significant share in the sequence.

This observation proves that the attention mechanism can successfully capture
the information that is helpful to the downstream tasks in the aggregation of
paths.

##  6 Conclusion

In this paper, we propose a model—PathSAGE, which can expand the size of
receptive field without stacking layers which solved the problem of “neighbor
explosion” and “over smoothing”. Regarding all paths as the neighbor nodes of
the central node, PathSAGE samples paths starting from the central node based
on random walk and aggregates these features by a strong encoder—Transformer.
Consequently, the model can obtain more features of neighbor nodes and more
patterns around the central node. In our experiment, PathSAGE achieves the
state-of-the-art on inductive learning tasks. Our model provides a novel idea
to handle large-scale graph data.

## Acknowledgments

. This work is supported by National Key R&D Program of China
(No.2017YFB1402405-5), and the Fundamental Research Funds for the Central
Universities (No.2020CDCGJSJ0042). The authors thank all anonymous reviewers
for their constructive comments.

## References

  * [1] Bruna, J., Zaremba, W., Szlam, A., LeCun, Y.: Spectral networks and locally connected networks on graphs (2014)
  * [2] Chen, J., Zhu, J., Song, L.: Stochastic training of graph convolutional networks with variance reduction (2018)
  * [3] Chen, J., Ma, T., Xiao, C.: Fastgcn: Fast learning with graph convolutional networks via importance sampling (2018)
  * [4] Chiang, W.L., Liu, X., Si, S., Li, Y., Bengio, S., Hsieh, C.J.: Cluster-gcn: An efficient algorithm for training deep and large graph convolutional networks. In: Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining. pp. 257–266 (2019)
  * [5] Defferrard, M., Bresson, X., Vandergheynst, P.: Convolutional neural networks on graphs with fast localized spectral filtering (2017)
  * [6] Gao, H., Wang, Z., Ji, S.: Large-scale learnable graph convolutional networks. In: Proceedings of the 24th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining. pp. 1416–1424 (2018)
  * [7] Hamilton, W.L., Ying, R., Leskovec, J.: Inductive representation learning on large graphs (2018)
  * [8] Hochreiter, S., Schmidhuber, J.: Long short-term memory. Neural computation 9(8), 1735–1780 (1997)
  * [9] Huang, W., Zhang, T., Rong, Y., Huang, J.: Adaptive sampling towards fast graph representation learning. arXiv preprint arXiv:1809.05343 (2018)
  * [10] Kingma, D.P., Ba, J.: Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6980 (2014)
  * [11] Kipf, T.N., Welling, M.: Semi-supervised classification with graph convolutional networks (2017)
  * [12] Klicpera, J., Bojchevski, A., Günnemann, S.: Predict then propagate: Graph neural networks meet personalized pagerank. arXiv preprint arXiv:1810.05997 (2018)
  * [13] Li, Q., Han, Z., Wu, X.M.: Deeper insights into graph convolutional networks for semi-supervised learning. In: Proceedings of the AAAI Conference on Artificial Intelligence. vol. 32 (2018)
  * [14] Liu, Z., Chen, C., Li, L., Zhou, J., Li, X., Song, L., Qi, Y.: Geniepath: Graph neural networks with adaptive receptive paths. In: Proceedings of the AAAI Conference on Artificial Intelligence. vol. 33, pp. 4424–4431 (2019)
  * [15] Monti, F., Boscaini, D., Masci, J., Rodola, E., Svoboda, J., Bronstein, M.M.: Geometric deep learning on graphs and manifolds using mixture model cnns. In: Proceedings of the IEEE conference on computer vision and pattern recognition. pp. 5115–5124 (2017)
  * [16] Paszke, A., Gross, S., Chintala, S., Chanan, G., Yang, E., DeVito, Z., Lin, Z., Desmaison, A., Antiga, L., Lerer, A.: Automatic differentiation in pytorch (2017)
  * [17] Vaswani, A., Shazeer, N., Parmar, N., Uszkoreit, J., Jones, L., Gomez, A.N., Kaiser, L., Polosukhin, I.: Attention is all you need (2017)
  * [18] Veličković, P., Cucurull, G., Casanova, A., Romero, A., Lio, P., Bengio, Y.: Graph attention networks. arXiv preprint arXiv:1710.10903 (2017)
  * [19] Yang, Y., Wang, X., Song, M., Yuan, J., Tao, D.: Spagan: Shortest path graph attention network (2021)
  * [20] Ying, R., You, J., Morris, C., Ren, X., Hamilton, W.L., Leskovec, J.: Hierarchical graph representation learning with differentiable pooling. arXiv preprint arXiv:1806.08804 (2018)
  * [21] Zeng, H., Zhou, H., Srivastava, A., Kannan, R., Prasanna, V.: Accurate, efficient and scalable graph embedding. In: 2019 IEEE International Parallel and Distributed Processing Symposium (IPDPS). pp. 462–471. IEEE (2019)
  * [22] Zeng, H., Zhou, H., Srivastava, A., Kannan, R., Prasanna, V.: Graphsaint: Graph sampling based inductive learning method (2020)
  * [23] Zhang, J., Shi, X., Xie, J., Ma, H., King, I., Yeung, D.Y.: Gaan: Gated attention networks for learning on large and spatiotemporal graphs (2018)
  * [24] Zhao, Z., Chen, H., Zhang, J., Zhao, X., Liu, T., Lu, W., Chen, X., Deng, H., Ju, Q., Du, X.: Uer: An open-source toolkit for pre-training models. EMNLP-IJCNLP 2019 p. 241 (2019)
  * [25] Zhou, J., Cui, G., Zhang, Z., Yang, C., Liu, Z., Wang, L., Li, C., Sun, M.: Graph neural networks: A review of methods and applications. arXiv preprint arXiv:1812.08434 (2018)