Stefan Problem through Extended Finite Elements: Review and Further Investigations

A general review of the extended finite element method and its application to the simulation of first-order phase transitions is provided. Detailed numerical investigations are then performed by focusing on the one-dimensional case and studying: (i) spatial and temporal discretisations, (ii) different numerical techniques for the interface-condition enforcement, and (iii) different treatments for the blending elements. An embedded-discontinuity finite element approach is also developed and compared with the extended finite element method, so that a clearer insight of the latter can be given. Numerical examples for melting/solidification in planar, cylindrical, and spherical symmetry are presented and the results are compared with analytical solutions.


Introduction
In its original formulation, the Stefan problem is a representation of (temperature-driven) first-order phase transitions in matter through a boundary value problem for PDEs in which a discontinuity surface, internal to the domain, can move with time.Indeed, the problem was originally formulated in relation to ice formation and melting.However, the same mathematical framework can be applied to a variety of problems arising in different areas of science and technology, such as the diffusion of gasses in biologic tissues, biofilm and hydrogel growth [1,2], the penetration of solvents in polymers, the OPEN ACCESS flow in porous media, filtration problems, free surface flows, etching, shock propagation, and financial mathematics (e.g.[3]).For the formulation details of some of the physical problems mentioned above, and for further examples, see Reference [4].
Since the topic of this paper is the numerical treatment of the problem, we select a specific physical case, that is the solid-liquid phase transition (e.g.[5]).Applicative examples are the casting of metals, the freezing and thawing of foods, the production of ice, or ice formation on pipes.As in all the firstorder phase transitions, a discontinuity in the first derivative of the free energy with respect to a thermodynamic variable is exhibited and a latent heat is involved.In a solid-liquid transition at constant volume, the pressure, i.e. the first derivative of the free energy with respect to the specific volume, is discontinuous.During the transition, the system either absorbs or releases a fixed amount of energy at a constant temperature.Since the heat cannot be exchanged instantly, a front of solidification or melting is present and moves according to the heat transport across the front itself.
From the mathematical point of view, the solid and liquid regions are subdomains in which the coefficients (representing the physical properties of the medium in each phase) of the underlying PDEs are continuous and differentiable up to the order of the PDEs.The coefficients are discontinuous across the surfaces that separate the adjacent phases and the PDEs are not valid there, so that additional equations are needed for closure.These are derived from energy conservation, namely by the Stefan condition that expresses the local velocity of the moving freezing/melting front as a function of the heat flux evaluated at both sides of the phase boundary.Therefore, the solution of the heatconduction equations are required within unknown subdomains (the liquid and the solid regions) and their interface must be determined as part of the solution.
Existence and uniqueness of the solution to Stefan problem are demonstrated in References [6] and [7], respectively.However, analytical solutions are available in a close form only for a restricted number of simple, particular cases, all characterised by a high degree of symmetry in the geometry and in the boundary and initial conditions.For all the other cases a numerical treatment is required.
For developing a numerical scheme for Stefan problem, two issues must be solved.Firstly, the evolving geometry of the phase-interface must be suitably described with a discrete model; secondly, the temperature field must be discretised in such a way as the jump in the heat flux, i.e. in the temperature gradient, at the phase interface can be reproduced.The description of the phase-interface can be done explicitly or implicitly.
In the explicit methods, a Lagrangian approach is followed and some marker points (finite in number) on the phase interface are selected and tracked explicitly.This is the case, for instance, in which the temperature field is discretised by standard Finite Element Method (FEM); the mesh is generated at the initial time instant in such a way as the element boundaries lay on the phase interface.The finite element nodes on the interface are used as marker points and, at each time step, are displaced according to the solution of the heat-transfer problem within each phase.Eventually, the initial FE mesh will become too distorted and a new mesh must be generated in the solid and liquid regions.
In the implicit methods, on the other hand, the computational grid is usually fixed and the position of the interface is obtained indirectly from some field defined over the whole domain (and suitably discretised).Examples of implicit methods are the enthalpy method, the phase field method, and the level set method.In the enthalpy method (e.g.[8,9]), the enthalpy field is considered and the position of the interface is determined in a smeared way by a jump of such field.In the phase field method (e.g.[10]) a phase function is defined over the domain; it assumes fixed values in each of the phases (e.g. 1 in solid and -1 in liquid) and varies smoothly between these values in the interface region.In order to allow such smooth variation, an artificial 'interface thickness' must be added to the model.Hence, the interface position is not defined exactly, although it can be chosen conventionally as the surface where the phase field assumes an average value (e.g.zero) or as some kind of 'average surface' of the interface region; hence, a sufficiently fine mesh is required in order to resolve the interface zone.Of course, it is necessary to provide additional governing equation for the phase field evolution in a thermodynamically-consistent way.Finally, in the level set method a function is defined all over the domain and the position of the phase interface is obtained as a level set of this function.In comparing the last two methods, we see that, although the level set function is artificial, the physics of the problem is fully respected; on the other hand, the phase field has a physical meaning, but it necessitates of an artificial 'thickness' of the interface.Some comparative analyses in Reference [11] lead to the conclusion that methods based on the level set are the most promising and general.For these reasons, and since it allows a general, elegant formulation of the extended finite element method, the level set method is chosen here for the numerical analyses and it is reviewed in Section 3.
The second issue in the numerical treatment of the Stefan problem is the discretisation of the temperature field.The crucial feature, in order to account for the evolution of the phase interface, is the discontinuity in the temperature gradient across the interface itself.
In the aforementioned enthalpy method the problem is reformulated in terms of enthalpy instead of temperature.Although this approach has the advantage of offering a straightforward discretisation of the problem, it has two main disadvantages, as the discrete set of equation in terms of enthalpy becomes nonlinear and the phase interface is smeared throughout a finite thickness.In the so-called boundary immobilisation method, a new set of coordinates is defined in which the phase-boundary is fixed.Caldwell and Kwan, in their comparative analysis [12], highlight the effectiveness of this method for one-dimensional problems, but discard it for planar or spatial cases because of the difficulties in defining the coordinate transform with respect to a reference coordinate system fixed in time.Other methods, which are not discussed in depth here due to their lack of generality, are the perturbation method, which requires hints from a partial analytical solution, the nodal integral method, which only applies to plane symmetric geometries, and the heat balance integral method, which is restricted to time-independent boundary conditions.For details on these methods, we remand to the last quoted reference and further references therein.
More general approaches to the temperature field discretisation are based on standard techniques, such as the finite difference method with a moving grid (e.g.[11]), the boundary element method, and, above all, the finite element method with adaptive mesh (for references on the last two methods applied to Stefan problem, see in Reference [13]).As already mentioned while discussing the interface tracking methods, Stefan problem can be treated by using standard FEM and a mesh which conforms to the phase interface (e.g.[14]).The heat-conduction problem is then solved separately in the solid and liquid domains by imposing essential boundary conditions on the phase-interface, i.e. by prescribing the phase-transition temperature.Then, the evolution of the interface is obtained from the jump in the heat flow across the interface, and the position of the nodes belonging to the interface is updated accordingly.The jump in the temperature gradient across the phase interface is possible thanks to the adapted mesh but, of course, remeshing is required every few discrete time-steps.Moreover, each time a new mesh is generated, the nodal values of the temperature must be recomputed by using the interpolated values from the older mesh.
In order to avoid remeshing it is possible to enhance the FEM by allowing for discontinuities that cross the elements.This can be achieved either by embedding the discontinuity directly within the finite elements intersected by the phase-interface and controlling it by additional elemental degrees of freedom, or by enriching the interpolation by means of nodal extended shape functions controlled by additional nodal degrees of freedom.The former approach we call embedded-discontinuity FEM, while the latter is well known as eXtended Finite Element Method (XFEM).In this paper, we focus on these two methods by reviewing the XFEM and adapting it to a specific case, and by developing an original embedded-discontinuity FEM.In both methods, the nodes of the discretisation are fixed and no remeshing is required.It is then natural to adopt an Eulerian approach for the interface description, so that the level set method is used.
The XFEM has been applied to Stefan problem by Merle and Dolbow (2002) [15], who focus on the one-dimensional case, by Chessa et al. (2002) [16], who develop the general method that is adapted here for the use with different constraint enforcement methods, and by Zabaras et al. (2006) [17], who introduce a simplified case of fluid motion and dendritic solidification by using a smeared-interface model.An application of XFEM to the problem of biofilm growth in two dimensions is developed in Reference [18].No applications of the embedded-discontinuity FEM to Stefan problem are known to the writers.
Though, as it has been shown, the research is currently tackling planar and spatial cases, the onedimensional Stefan problem is the subject of several recent research works [11,12,15] and of a monographic book [19].Also the examples that are presented in this paper are developed for cases that can be treated as one dimensional due to planar, cylindrical, or spherical symmetry.This is done for the sake of simplicity in exposition, implementation, and result visualisation, which allow concentrating on the differences between different approaches of the XFEM and between XFEM and embedded-discontinuity FEM.Moreover, exact analytical solutions are available in some onedimensional cases, so that error and convergence can be studied.In further cases, it is possible to develop, as it is done in Sections 6 and 7, simple approximated analytical solutions in the so-called pseudo steady state approximation, in which the heat-conduction problem and the evolution of the interface can be decoupled.Though the aim of the present one-dimensional analyses is basically heuristic, they may also find practical applications.Examples in planar, cylindrical, and spherical symmetry are, respectively, the maintenance of ice surface of skating facilities, the freezing of water into pipelines, or the solidification of droplets of liquid metal into cooling fluid.The initial review provided in this paper includes details that are not strictly necessary for the subsequent onedimensional applications but that are advisable for a clearer understanding and propaedeutical for the discussion.
This paper is organised as follows.In Section 2, the governing equations of the melting/solidification problem are briefly recalled.Sections 3 and 4 are review sections, dedicated to the level set method for the geometrical description of the phase interface and its evolution, and to the extended finite element method, respectively.In Section 5, a method for embedding the gradient discontinuity into finite elements is developed and a specific element is formulated.In Section 6, the XFEM is thoroughly discussed by means of numerical applications to planar-symmetric problems; spatial and temporal discretisation errors, the enforcement of the interface condition, and the treatment of the so-called blending elements are studied; moreover, the XFEM is compared to the embeddeddiscontinuity FEM.In Section 7, polar-symmetric problems are analysed and the relevant XFEM results are compared to analytical solutions.

Governing Equations
The Stefan problem of melting/solidification in absence of material transport is considered.For the case of the melting/solidification of a two-phase flow, see Reference [20].
Let  be the domain shown in Figure 1 and  its boundary with outward normal unit vector n .The domain is composed of a solid region sol  and a liquid region liq sol \     .We denote by tran  the interface between the two regions and by sol n the outward normal from sol  .The temperature in The model accounts for the heat transferred by conduction only; therefore, convection, thermal expansion, and buoyancy are not considered.Since deformations are excluded, the mass density  is the same in the solid and liquid phases.The energy conservation equation in the domain reads where c is the heat capacity, q the heat flux vector, and s the volumetric heat source.
The constitutive equation which characterises the conductivity within the material is given by Fourier law, where k is the thermal conductivity second-order tensor.In the following, we assume isotropy with respect to heat conduction, so that k  k 1, where 1 denotes the identity tensor and k the scalar thermal conductivity.k and c are taken as constants in each phase and their values in the solid and in the liquid are indicated by subscripts 'sol' and 'liq' respectively.
The essential and natural boundary conditions are ˆon , ˆon , where   ˆ, T t x and   ˆ, q t x are assigned functions defined over respectively, with T q     and T q     .
The phase-interface is characterised by the constant temperature tran T , so that the condition holds.It is noted that the temperature field must be continuous through the interface; this requirement must be included (or enforced) in the numerical models as discussed in Section 4.
It is also possible to consider more complicate interface conditions which involve the velocity of the interface or its geometrical properties.An example is Gibbs-Thomson relation, which describes the unstable dendritic growth of crystals into an undercooled melt, and modifies Eq. ( 4) as , where   and v  are the surface tension and the kinetic mobility coefficients, sol div   n is the mean curvature of the interface, and v its normal speed (negative for melting and positive for solidification).
The evolution of the interface is governed by the Stefan equation, which expresses the energy balance of the interface, where L is the latent heat of the phase transition and   q the jump in the heat flux normal to tran  given by The initial condition for the transient problem is provided by the knowledge of the temperature field at time 0 where 0 T is a known continuous function over  .The initial conditions are completed by the knowledge of the position of the interface at 0 t  .

Description of the interface through level set
The level set method (e.g.[21]) is a numerical scheme used to describe implicitly geometrical objects such as surfaces and lines.The basic concept is to represent a surface  , of co-dimension one, in the domain  as the isocontour of some function    x defined over  .The method can be extended to model objects with higher co-dimensions by considering multiple level set functions (as many as the co-dimension of the object); for instance, a curve in the three-dimensional space can be represented as the intersection of two surfaces, i.e. the isocontours of two functions   where the arbitrary sign convention in which  is positive in the liquid, negative in the solid, and vanishes on the phase interface, is adopted.Eq. ( 9) allows the construction of  , when the position of the interface is known.Vice versa, when  is given, the interface tran  is the zero isocontour of  , Take note that, whereas the function   , d t x is not differentiable in space across the interface, the gradient of   ,t  x is defined everywhere in  (provided that the interface is sufficiently regular) and on the interface it points in the direction of the outward normal from the solid region.Therefore, although the interface is obviously the zero level set of d , d itself is not suitable for describing the interface position since it is not sufficiently regular and, in the numerically approximated form used in computational models, would not reach the zero level in certain parts or would cross the zero twice in others.
Besides its position, the level set allows one to describe further geometrical properties of the interface.The normal unit vector is given by sol Since the level set function has been chosen as the signed distance function, 1    , so that Eq. .The interface mean curvature, which may have a physical interest in the case of Gibbs-Thomson diffusion (see Section 2), can be computed as sol div div which simplifies to 2     in the case of signed distance function.
In order to track the movement of the interface it is necessary to define suitable evolution equations for the level-set update.An Eulerian approach is followed by letting the level-set field be advected by a suitable velocity field v .The relevant evolution equation is then given in the Hamilton-Jacobi form by requiring the material derivative of the level set to vanish, In particular, Eq. ( 13) requires that  be constant on the interface.The field   ,t v x represents an extension of the interface velocity, which is physically defined only on tran  , to the whole domain  .
In fact, whereas in a Lagrangian approach, where the position of the interface is tracked explicitly, the velocity on tran  suffices for its evolution, in an Eulerian approach, it is necessary to extend the velocity field to the entire domain so that the level set function, which controls the interface position and is defined everywhere in  , can be advected by such extended velocity field.It is convenient to introduce the scalar normal velocity that is the component of v in the direction of   ; Eq. ( 13) rewrites then as The field F is the extension of the normal interface velocity to the entire domain, and it is known only on the interface, where it is physically meaningful and can be obtained by Eq. ( 5), Since the velocity field in tran \   is artificial, there is a certain freedom in its construction, the only requirements being that F be continuous through the interface and satisfy Eq. ( 16).It can be demonstrated [22] that, if F and  are smooth and their gradients are orthogonal, then F tends to preserve the initial signed-distance properties of the level set.
The evolution of   ,t  x is controlled by Eq. ( 15), with the initial condition Practically, in the initial time instant, 0  is constructed by the knowledge of   tran 0  by using the definition in Eq. ( 9), whereas in the subsequent time instants, the updated position of the interface is obtained from   t  by Eq. (10).The field F is constructed by solving Eq. ( 17) with the essential boundary conditions provided by Eq. (16).The problem can be casted in a discrete form by using standard finite element techniques developed from the weak forms of Eqs ( 15) and ( 17), where  and F  are test fields belonging to   1 H  .Essential boundary conditions are provided by a discretised version of Eq. ( 16), which can be obtained, for instance, by considering the projection of the interface velocity onto the finite element nodes whose support is intersected by the interface.The complete procedure is thoroughfully derived in Reference [16], where the numerical issues arising in the integration of first order hyperbolic equations are solved by adding stabilisation and shockcapturing terms; further stabilisation techniques for the update of the level set function are described in Reference [21].In the numerical analyses, the discretised field  progressively diverges from a signed distance function, so that every few steps a reinitialisation must be performed; some relevant procedures based on the fast marching algorithm are described in References [22] and [17].

XFEM approximation
The basic concept within the XFEM is to locally enrich the finite element approximation so that certain features of the considered problem can be reproduced, which are precluded to the standard FEM (e.g.[23]) unless ad-hoc elements or mesh adaptivity are used.    for 1 2 e e  ; and a finite number n of points, the FE nodes, whose position is denoted by is used for controlling the temperature field.The standard finite element approximation of the temperature field is where i N are the standard finite element shape functions and the parameter i T represents the value of the temperature at node i .In fact, the basis functions i N are constructed in such a way that , being ij  Kronecker delta, and Figure 2 shows an example of 2D shape functions and highlight the support of the i-th shape function, denoted as supp i N .If Eq. ( 21) holds, the set of functions   i N is called a partition of unity of  .
Additional extended basis functions j  , 1, , j n    , can be tailored on the specific problem so that they are able to reproduce the desired local feature, e.g. a discontinuity.In order to recover conformity, i.e. continuity between neighbouring elements, the extended basis functions are multiplied by a partition of unity, which in FE-based approaches is naturally provided by the set of the FE shape functions, see Eq. (21).A general form of the enriched interpolation is where j i a are additional nodal degrees of freedom.Eq. ( 22) characterises the Partition of Unity Finite Element Method (PUFEM, e.g.[24]) also called Generalised Finite Element Method (GFEM, e.g.[25,26]), which combines the FEM and the Partition of Unity Method (PUM), introduced in References [27,28].For a recent review see also Reference [29].The PUM has been also combined with meshless approaches (e.g.[30][31][32][33]), resulting in enriched meshfree methods (mostly applied to fracture mechanics, e.g.[34][35][36][37][38]), and with smoothed finite elements [39], which are less sensitive than standard finite elements to mesh distortion and represent an interesting approach in case of remeshing (e.g.[40][41][42]).
The extended finite element method (e.g.[43][44][45][46][47][48][49]) can be regarded to as a local version of the PUFEM.In fact, the special feature to be capture, namely a discontinuity, is usually localised; therefore, the extended basis is useful only in a restricted region surrounding the discontinuity itself.The subset of shape functions, whose support includes one specific element, is a partition of unity on that element.Therefore, it is sufficient to consider the shape functions relevant to the nodes whose support is intersected by the discontinuity  and Eq. ( 22) can be rewritten as where We define the union of the element intersected by the discontinuity as and the union of the support of the enriched nodes as supp with In order the nodal parameters i T to maintain the meaning of nodal temperatures they have in standard FEM, Eq. ( 23) can be further rewritten as where It is straightforward to observe that   , 0 holds in Eq. ( 29) as in the standard FE approximation of Eq. (20).Although the original enrichment given by   1 , ( ) is modified and does not simply superimpose to the standard FE approximation as in Eq. ( 23), its most interesting property, i.e. the ability of representing discontinuity, is kept.
The type and number of enrichment functions must be chosen according to the physics being considered: Heaviside function (or some regularised version of it) is the natural choice in case of strong discontinuities across an interface such as a crack; suitable enrichment functions can be defined around crack tips too (e.g.[43]).In the Stefan problem a discontinuity of the temperature gradient normal to the interface must be modelled.One of the simplest functions reproducing such gradient discontinuity (and preserving continuity) is the absolute value function.
The extended finite element method is easily coupled with the level set method, as the latter provides a practical and general way for representing the location of the discontinuities.For this reason the extended basis functions are expressed in terms of the level set function.In the case of Stefan problem,  is the phase-transition interface tran  and the gradient discontinuity in the direction normal to the interface is represented by using 1 Then, Eq. ( 28) simplifies in the form with In this way, the temperature field is approximated by using as basis functions the union of the standard finite element shape functions with a set of enrichment functions.Take note that since the position of the interface evolves with time, the set of enrichment function is time-dependent too.
Eq. ( 31) can be rewritten in matrix form as where the matrix collects the standard (timeindependent) and the enriched (time-dependent) shape functions, and the vector , the standard and additional nodal degrees of freedom.For all that is discussed above, this choice of enrichment function, initially proposed by Chessa et al. [16], seems the most consistent with the viewpoint of the XFEM.
A different XFEM model for Stefan problem is proposed for the one-dimensional case in Reference [15], where Heaviside function is chosen for the enrichment and the device in Eq. ( 29) is not used, so that . This approach does not guarantee a priori the continuity of the temperature field across the interface, which is then (weakly) enforced.In any case, the continuity can be restored exactly only in one-dimensional cases.The choice of Heaviside function, which is natural in the case of strong discontinuities (we remind that the XFEM has been initially and intensively developed for fracture mechanics, see References [50,51] and the review in Reference [52]), seems complicated in the case of the gradient discontinuity required by Stefan problem.Moreover, with those enrichment functions, the parameters i T lose the meaning of nodal temperatures, since one has that , which is in general different from   i T t .This issue could be cured by using the idea in Eq. ( 29) and rewriting the enrichment as Operatively, in the XFEM the usual FE mesh is firstly produced and then, by considering the location of the discontinuities, a few nodal degrees of freedom are added.Only the nodes whose support is intersected by the interface are enriched, as shown in Figure 3.For implementing this condition, Zabaras et al. [17] suggest to enrich a node if at least one of the element edges containing it is intersected by the interface (this criterion hold of course only in the case of finite elements with no internal nodes).Since the interface evolves with time, the set of enriched elements and nodes is also modified during a numerical simulation.At a given time instant, it is possible to distinguish three types of elements, as depicted in Figure 3: (i) elements with all nodes enriched, (ii) elements with only some nodes enriched, and (iii) elements with no enriched nodes.The first group is formed by the elements intersected by the phase-transition interface and their union is   .In the region   the enrichment is fully considered and conformity is ensured, as the enrichment is multiplied by a partition of unity, see Eq. ( 27).The union of the elements in the second group is \ A    .The treatment of these elements is discussed in detail in Section 6.5; here, it is noted that they do not have full enrichment since the subset A of the standard shape functions does not constitute a full partition of unity in this region.These elements are often referred to as blending elements, since, there, the interpolation 'blends' from the enriched to the non-enriched one.
As to the elements of the third group, their union is \ A   .In this region, according to Eq. ( 28), the enrichment is not accounted at all and only the standard FE approximation holds.

Semi-discrete equations
In order to produce a discrete model, besides introducing an approximation as discussed in the previous Section, it is necessary to weaken the continuum problem.Here, the integral version of the energy balance is used.The weak form of Eq. ( 1), with the boundary conditions of Eq. ( 3) and the substitution of the isotropic version of the constitutive model in Eq. ( 2), reads Eq. ( 34) must be discretised in space, through Eq. ( 33), and time.Since the shape function matrix N depends on time through the position of the interface, the time derivation of Eq. ( 33) would require deriving the extended shape functions and, in turn, the level set.It is therefore more convenient to discretise first in time and then to substitute the spatial discretisation.By using backward Euler scheme (e.g.[53]), the semi-discrete version of Eq. ( 34) is obtained, in which the superscript notation in parentheses denotes the discrete time steps, i.e.
, p   , and t  is the selected time step (the symbol  is used here to denote a finite difference, as usual in Computational Mechanics).Finally, the discrete form of the problem is obtained by substituting Eq. ( 33) into Eq.( 35), where For the sake of brevity, we rewrite Eq. ( 36) as with By solving the linear system of algebraic equations in Eq. (38), it is possible to obtain the discrete temperature configuration   p T at the current time step, given the one   1 p

T
at the previous time step.In order to evaluate   p M  it is necessary to know also the level set at both the current and the previous time step.The overall procedure can be explicit, if the level set update is also explicit, i.e. if the level set at current time step can be obtained from the temperature field and level set at the previous time step (see Section 3 and References [16,17]).
On the other hand, if an implicit time stepping is used, an iterative procedure is required for the simultaneous solution of Eq. ( 38) and the update of the level set.This technique is used by Merle and Dobow [15], who use a generalised one-step trapezoidal integration algorithm with Finally, we would like to remark that the evaluation of the matrixes and vector in Eq. (37) requires a careful numerical integration, which must take into account the presence of the discontinuity in the enriched elements.Several algorithms have been developed in the literature (see Reference [54] and further references therein).Examples are Gaussian integration into the sub-elements obtained from the enriched elements cut by the interface, and the integration by subdivision into simplexes (triangles in 2D or tetrahedra in 3D).In the analysis below, the algorithm suggested in Reference [17], is used, which is based on rectangular rule and regular grids of sampling points.This integration technique is simple to implement, though relatively expensive from the computational point of view.

Phase-interface temperature condition
Eq. ( 38) represents the discrete form of the problem of heat conduction as from Eqs (1-3).The interface temperature condition expressed by Eq. ( 4) must be also recasted into a discrete form and enforced.This can be achieved by following two different approaches: (a) a general one in which the enforcement methods (see below) are applied directly to Eq. ( 4), a weak form is derived, and, then, discretised; and (b) a simplified approach in which Eq. ( 4) is firstly discretised and then the constraint enforcement methods are applied to the discrete form.The simplified approach (b) is followed here.This can be successfully used for the one-dimensional cases presented in the examples of this paper.However, the enforcement of interface conditions in methods derived from the PUM, such as the XFEM, is an intensively-researched field, due to the stability issues arising in higher dimensional cases (e.g.[55,56]).
Eq. ( 4) is discretised by selecting some points, finite in number, on the interface, namely   , at time step p .By using Eq. ( 33) one obtains The points where the interface condition is imposed can be individuated as the intersections of the interface with the edges of the finite element mesh.At a given time step, Eq. ( 40) represents a linear constraint, which can be rewritten in matrix form as where the all components of the vector c are equal to tran T and the h-th row of the matrix , which can be determined from the knowledge of the level set.Two general numerical methods for enforcing the constraints are considered here: (i) the penalty method and (ii) the Lagrange multiplier method.Although these methods can also be applied to nonconservative problems [57], they are presented here in their conservative version.
The idea within penalty method is to add to the potential, of which Eq. ( 38) represent the stationary condition, an artificial term 1 2 T bg g , where b is a positive constant, in general dimensional, called penalty parameter.If b is sufficiently large, the minimum of the modified potential cannot be achieved without g being very small, i.e. without (approximatively) satisfying the constraints.The stationary condition of the additional potential provide the penalty term

 
T b  G GT c to be added to the righthand side of Eq. (38), so that the temperature update step becomes where The solution generated with the penalty method satisfies the condition over the interface only approximately.The choice of the parameter b is crucial, since a value 'too small' will lead to unacceptable looseness of the constraints, whereas a value 'too large' will deteriorate the conditioning of the linear problem.Examples of this are provided in Section 6. 4, where a suitable dimensionless version of b is defined in order to perform problem-independent analyses.
In the Lagrange multiplier method, the constraints are appended to the potential minimisation problem through Lagrange multipliers, which become additional unknowns.A general approach and a stabilisation proposal are discussed in Reference [58].In the simplified approach followed here, Lagrange multipliers are applied to the discretised constraint of Eq. ( 40), so that that they belong to their dual finite-dimensional vector-space.The term added to the potential is T g l , where l is a vector collecting the discrete Lagrange multipliers, and the term T G l must be added to the right-hand side of Eq. (38).The potential stationary condition is sought also with respect to l , so that additional equations, i.e. the constraints themselves, become part of the system, The disadvantages of the Lagrange multiplier method are the increased dimension of the system and its possible loss of positive-definiteness.

Finite elements with embedded gradient discontinuity
An alternative approach to model localised discontinuities within the FEM consists in embedding the discontinuity directly into the elements by adding degrees of freedom at elemental level instead of nodal level (as in the XFEM).The resulting methods appear in the literature with different names and they can be considered somehow as precursors of the XFEM (e.g.[59,60]).Further references, mainly focusing on fracture mechanics, can be retrieved in the comparative study by Jirásek [61].
The basic idea is to use standard FE interpolation as from Eq. ( 20) in all the elements but those intersected by  .In these latter, the standard interpolation is enriched by suitable functions ˆi ), which vanish at the nodes and reproduce the desired discontinuity.The resulting interpolation can be written as where ˆi a are element-level parameters.The enrichment functions ˆi  depend on the time through the interface-position, which, in turn, is described by the time-evolving level set function.Eq. ( 45) is the embedded-discontinuity version of XFEM's Eq. ( 31) and can be rewritten in matrix form as where the matrix collects the standard and discontinuity-embedding shape functions, and the vector the standard nodal degrees of freedom and the additional elemental ones.
The numerical model is developed by following the same steps as in Sections 4.2 and 4.3, so that the same formal expressions obtained there for the XFEM hold for the embedded-discontinuity FEM, once Eq. ( 31) is replaced by Eq. (45).
In order to explain the procedure, a two-node element with an embedded gradient discontinuity for the one-dimensional Stefan problem (see also Section 6) is developed here.In this simple case the standard FE shape functions are linear (Figure 4a) and the enrichment of the element containing the phase interface can be obtained by a piecewise-linear hat-shaped function   ˆ, x t  (Figure 4c), whose contribution vanishes at the nodes, assumes a unit value at the discontinuity, and it is controlled by the elemental degree of freedom â .It can be observed that the definition of  given above does not cover the case in which the phase interface coincides with a node.When the interface reaches the local i-th node ( 1, 2 i  ), the function  coincides with i N , in the element that last contained the interface.However,  is still linearly independent from i N at the global level, as visualised in Figure 5, and no particular indeterminacy arises, unless the interface reaches the domain boundary.In the last, quite particular case, the loss of linear independency may also compromise the convergence of the gradient discontinuity at the phase boundary, as observed by Merle and Dolbow [15], and, in any case, the conditioning of the system matrixes deteriorates as the phase interface gets closer to a boundary node.The case discussed here is particularly simple and further difficulties arise in two and three dimensions or in case of stronger discontinuities, as in fracture mechanics.For these reasons, XFEM is usually preferred to embedded-discontinuity methods in those fields.However, in the case of onedimensional Stefan problem, the approach developed above proves quite effective, as shown by the examples in Section 6.6.

Numerical study on the error and applications to planar-symmetric problems
Exact solutions to Stefan problem exist in a very small number of cases such as those of planar symmetry with semi-infinite domain, in which the interface position is proportional to the square root of the time (similarity solutions).These cases are reviewed in this Section, so that the numerical results can be compared to exact values.The accuracy of XFEM is evaluated with respect to space and time discretisations.Penalty and Lagrange multiplier methods to enforce the interface temperature condition are compared, and the treatment of the blending elements is discussed.XFEM and embedded-discontinuity FEM are compared too.Finally, a second planar-symmetric example with finite domain is presented.
Besides exact solutions and for the cases in which such solutions are not available, we also consider the so-called pseudo steady state (PSS) approximations.These are based on the assumption that the rate of the interface movement is much slower than the rate of the temperature diffusion within the domain, so that the temperature field at a given time instant can be approximated with the steady-state field corresponding to a fixed phase-interface, and the interface position can be, in turn, updated from the resulting steady-state temperature field.It follows that the PSS solution is a suitable approximation for large values of the Stefan number  , defined as where T  is some temperature difference, characteristic of the specific problem (e.g. the temperature difference between two thermostats providing essential boundary conditions).A PSS approximation provides a limit behaviour, useful for comparisons where no exact analytical solution is available, such as the cases considered in Sections 6.7 and 7.For further details on the exact and PSS analytical solutions, the Reader is referred to Reference [19].
From here on, planar-, cylindrical-, and spherical-symmetric problems are considered.These can be treated as spatially one-dimensional, allowing for some simplifications in the interface-position tracking.In particular, with the phase-transition interface reduced to one point in the one-dimensional space (representing a planar, cylindrical, or spherical surface in the original three-dimensional space), one coordinate suffices for individuating the interface position and generating the level set.
In the examples below, we consider one-dimensional extended finite elements with two-nodes and linear shape functions, enriched as in Eq. (31).The relevant enriched shape functions are shown in Figure 4b.In Section 6.6, we also consider the embedded-discontinuity FEM described in Section 5.Where it is not differently stated, the temperature condition on the phase interface is enforced through Lagrange multiplier method and no enrichment is used in the XFEM blending elements.

One-phase freezing problem
A planar-symmetric freezing problem in a semi-infinite domain is considered.Due to symmetry and material isotropy, the problem can be reduced to one dimension and the coordinate x in the direction orthogonal to the constant-temperature planes suffices to the problem description.The onedimensional domain With these initial and boundary conditions, no heat flow occurs in the liquid domain, which remains at a constant temperature, and the interface is driven only by the temperature gradient in the solid domain, so that the problem can be considered as a one-phase problem with moving boundary.
Neumann's exact solution is available for this case.The temperature field in the solid region is given by where the non-dimensional coefficient  is the positive real root of and the position of the freezing front is given by As further comparison, we also consider the PSS solution, although the Stefan number of the specific example, see Eq. ( 48), is not large enough for the PSS solution to constitute a good approximation [19].Under the PSS assumption, the one-dimensional governing equation becomes with the boundary conditions Integration of Eq. ( 52) leads to The further integration of Eq. ( 5), with and the initial condition Eqs ( 54) and ( 55) are the PSS versions of Eqs ( 49) and ( 51) respectively.All of them are written into a form in which mutual comparison is straightforward; in particular it is noted that the temperature profile becomes linear in the case of PSS solutions, whereas the expression of the position of the freezing front only differs on a time scale-factor.
In the numerical analyses, we use a finite domain . Since no heat flux occurs in the liquid phase as long as the position of the interface is contained into XFEM  , the numerical set-up is completely equivalent to the one studied analytically, provided that the natural boundary condition   ˆ0 q x    is imposed.In Figure 6, the interface position calculated through Neumann's and PSS analytical solutions are compared with the XFEM simulations obtained with a spatial discretisation in 10 equally-sized elements and a time-step 5 4 10 s t    . The numerical results are very close to the exact ones, whereas the PSS solution overestimates the freezing-front velocity because it does not account for the thermal inertia.At a closer look, the numerical solution lies slightly above the exact one, due to the use of backward Euler integration scheme, in which the interface position at current time-step is updated with the velocity of the previous time-step, i.e. a slightly overestimated one in this example.Of course, this effect is reduced with smaller time-steps and the XFEM solution converges to the exact one as the temporal and spatial discretisations are refined (see Sections 6.2 and 6.3).The temperature field at some different time-steps is shown in Figure 7, where the numerical and the exact solutions superimpose quite well one to the other, whereas the PSS approximation overestimates the phaseinterface positions.

Time discretisation
As mentioned above, backward Euler integration evaluates the interface position at every time step with the temperature gradient computed through finite differences in the previous time-step, so that the computed gradient is slightly larger, in the considered example, than the exact one.Figure 8 shows the interface position obtained with different time discretisations at a fixed time instant 6 10 s t  . The interface position converges to the exact one as t  goes to zero.The non-dimensional relative error is represented versus the time-step size in Figure 9.

Space discretisation
The dependence on the spatial discretisation of the error in XFEM simulations with respect to the exact solution is studied by considering the same example of Sections 6.1 with a fixed time step and a domain subdivision into different numbers of equally-sized elements.
The error in a spatial subdomain   which, for the convergence study, is equivalent to the square-root of the non-dimensional error in the energy.The convergence to exact solution is shown in Figure 10, where the error in the whole domain is averaged on the time and plotted as a function of the normalised element size (i.e. the inverse of the number of elements in which the spatial domain is subdivided).Of course, the error on the spatial discretisation is increased by the error on the interface position (and vice versa), which is responsible of the apparent scatter in the error values, as well as of a reduced convergence rate.In order to highlight this effect, the error on the spatial discretisation can be somehow 'depurated' from the timediscretisation error by enforcing the exact analytical interface position in the numerical simulations.In this way a higher convergence rate and no 'scatter' are obtained, as also shown in Figure 10.

Constraint enforcement
The enforcement of the phase-transition temperature at the interface, described in Section 4.3, is discussed here in the same settings of Section 6.1.
We consider first the penalty method, which is used in most of the research works on XFEM and Stefan problem [15][16][17], because of its robustness and simplicity of implementation.As discussed in Section 4.3, the main disadvantage of this method is that high values of the penalty coefficient b , while accurately enforcing the constraint, may lead to a poor conditioning of the numerical problem of Eq. (42).We define the error in the enforcement of the constraint at a given time instant as  cond [-] The conditioning of the linear system in Eq. ( 42) is measured by the condition number with respect to inversion, cond  , defined as the ratio between the maximum and the minimum singular values of the system matrix pen b  K K  ; the higher is cond  and the worse is the conditioning of the numerical problem.The error and the numerical conditioning are studied as functions of the non-dimensional penalty coefficient, which we define as with S and e  the average cross-sectional area and length of the one-dimensional elements, respectively.Since pen K is dimensionless, the term in parentheses in Eq. ( 59) contains the physical dimension and the average magnitude of the coefficients of the stiffness matrix K , which, in turn, contributes to K  .In this way,  is expected to be less problem-dependent than its dimensional version b .
Figure 11 shows the influence of the penalty coefficient on the quality of the constraint enforcement and of the numerical solution (both averaged over several time-steps); low values of  result in a poor enforcement of the interface condition, whereas high values lead to a progressive increase of the condition number.It is then evident that in practical cases the value of  must be chosen as a compromise between accuracy of the constraint enforcement and quality of the numerical solution.
The error in the constraint enforcement vanishes if Lagrange Multiplier method is used.The temperature field in the vicinity of the phase-interface at time 5 6 10 s t   as obtained with penalty ( 50   ) and Lagrange multiplier methods is shown in Figure 12 (position and temperature are indicated by solid dots for the nodes and by a circle for the interface).
In the analysed case, the condition number with the Lagrange multiplier method is 104.5 with an exact enforcement of the constraint; the same condition number can be obtained with the penalty method for 18   , with an error in the constrain enforcement exceeding 5%.Exact solution XFEM (Lagrange multiplier method) XFEM (penalty method =50)

On the interpolation in the blending elements
The treatment of blending elements (i.e.those whose nodal support is only partially enriched, see Section 4.1) is discussed here by means of a numerical example.Different approaches to the temperature interpolation in this region are possible: (i) the partial enrichment by routinely applying Eq. ( 31), (ii) the adoption of special strategies such as the use of hierarchical functions (e.g.[62]), and (iii) the use of standard FE shape functions only (as in the non-enriched elements).
The XFEM strategy consists in enriching the nodes whose support is intersected by the discontinuity.In this way the subset A of the standard FE shape functions constitutes a local partition of unity in the region   surrounding the discontinuity and does not affect the region far from the enriched nodes (where standard FE interpolation holds).However, by following approach (i), the enrichment involves a region A     and is only partial in the region \ A    , i.e. in the blending elements.In a blending element, only some nodes are enriched and the subset of the standard shape functions used to weight the discontinuity does not represent a partition of unity.This leads to the appearance of pathological terms and, in general, to a sub-optimal convergence rate [62].
One possible cure, the approach (ii), is the use of hierarchical, higher-order shape functions as proposed in Reference [62].The idea is that the enrichment functions are obtained by the product of the standard FE shape functions (linear, in our case) and the enrichment functions (also linear); therefore, it is possible to compensate their pathological contributions by using higher order FE shape functions in the blending elements (second grade polynomials).This approach proves very effective as demonstrated by their proposer.However, in those cases in which the interface is not fixed, local modifications of the topology are required as the interface evolves.When the interface moves from one element to another one, new nodes must be enriched and, consequently, some standard elements become blending elements.There, in a topological approach, new hierarchical nodes must be added.This may complicate the implementation of the method, even though the result is computationally cheaper than the remeshing required by standard FEM with adaptive mesh.Tarancón et al. [62], shows that the approximation error  in a two-node one-dimensional blending element satisfies the inequality where y and x are points inside the element and a is the extended degree of freedom of the enriched node.It is clear that the second term in the right-hand side of Eq. ( 60), which is not present in the standard FE case, may lead to an increased error.This pathological term can be compensated by using higher-order (second grade, in our case) standard shape functions for the blending element, controlled by additional standard degrees of freedom (hierarchical approach).Alternatively, as in the approach (iii), it is possible to eliminate the pathological term by considering the enrichment functions only in those elements whose support is completely enriched and simply using the standard FE shape functions in the blending elements.In fracture mechanics, where the XFEM has been originally developed, this approach causes problems, as the displacement field may not be well reproduced in the vicinity of the interface boundaries (the crack tips) and conformity between enriched and blending elements may be lost.In the case of Stephan problem, however, these are not concerns.In fact, the phase-interface has no boundaries inside the body (either it is a closed surface or its boundary belongs to the domain boundary) and the considered discontinuity is weaker (gradient discontinuity).Incidentally, we observe that non-conforming elements (sometimes called 'incompatible') are not uncommon, and the loss of conformity can be traded for an improved accuracy in secondary (but sometimes more relevant for technical applications) variables such as strains and stresses (e.g.[63][64][65]).
With that said, the approach (iii) seems optimal in the case of phase transitions and it is used in all the examples of the present paper.Only in this Section we compare the approaches (i) and (iii).Figure 13 shows the error as a function of the spatial discretisation, computed as in Section 6.3.The temperature fields obtained in the two cases at the same time instant are shown in Figure 14 (where the crosses denote the nodal positions and temperatures, and the circle the phase interface).It is noted that, not only approach (i) is less accurate than approach (iii), which confirms the estimates by Tarancón et al. [62], but it also leads to a poor estimate of the left-and right-gradients of the temperature at the interface, altering the rate of its movement.

Comparison with embedded-discontinuity finite elements
The one-phase problem of Section 6.1 is now analysed by means of the one-dimensional embeddeddiscontinuity FEM developed in Section 5.The simulations are performed by using the penalty method for the enforcement of the interface condition, with , the same discretisations used with XFEM in Section 6.1.Figure 6 shows the evolution of the melting front computed with embedded-discontinuity FEM in comparison with XFEM results and analytical calculations.The two numerical techniques provide comparable results, with only a slightly smaller accuracy in the embedded-discontinuity FEM.The temperature field in the element containing the discontinuity at 6 4 10 s t   is shown in Figure 15, where the two numerical methods appear as quite equivalent.A particular of the same plot is shown in Figure 16, where the differences in the interpolation (linear for embedded-discontinuity FEM and quadratic for XFEM) can be appreciated thanks to a strong magnification on the temperature axis.In Figure 15 and Figure 16, the dots denote nodal positions and temperatures and the circles highlight the phase interface.
As a second analysis, the convergence with spatial discretisation is studied in the same way of Section 6.3.The results are summarised by the plot in Figure 17, where embedded-discontinuity FEM and XFEM are compared.The two methods show, in the present case, comparable accuracy and convergence rate.

Two-phase solidification problem
The planar-symmetric solidification process of aluminium in a case in which the heat conduction occurs in both the solid and the liquid phases is considered.The one-dimensional domain is ; both phases are initially at the temperature 1 T .Under these conditions, no exact analytical solution is available.Then, Neumann's solution for the case in which the domain is unbounded on the side of positive x (i.e. when  tends to infinity) and the PSS solutions are used to compare the numerical results.Neumann's solution is exact in the semiinfinite domain; therefore, in the finite domain it provides a good approximation only in the earliest instants of the melting process.On the other hand, the PSS solution is expected to provide a better and better approximation as the system approaches the steady state.
For the two-phase problem described above, the interface position calculated through Neumann's analytical solution is where  gamma is the positive root of   the solid-to-fluid thermal diffusivity ratio.The PSS approximation is developed analogously to the one-phase case of Section 6.1.From Eq. ( 52) we obtain linear temperature fields in the solid and liquid phases, and the evolution of the melting front is obtained (in an implicit form) as where the constants      are introduced for the sake of conciseness.
The numerical data used in the example are reported in Table 1 and the resulting Stefan number is As it can be evinced from Figure 18, the interface positions evaluated by the XFEM and the two analytical models are close in the initial time instants.Afterwards, the numerical solution and the PSS approximation converge to the steady state, whereas Neumann's solution (which matches different boundary conditions) diverges.However, in the initial instants, when the position of thermostat 1 T makes little difference, Neumann's solution should be considered as reference, since it includes thermal inertia, contrarily to the PSS solution.The position of the phase-interface in such initial instants is represented in Figure 19, where the better agreement of the numerical simulation to Neumann's results is highlighted.
The temperature fields within the domain at some time instants are plotted in Figure 20.Initially, when the effects of thermal inertia are larger, the PSS (piecewise linear) temperature field poorly approximates the expected physical behaviour, which is more suitably represented by Neumann's solution.On the other hand, as the system approaches the steady state (represented by the dashed lines, almost superimposed with the fields at 20000 s t  ), the thermal inertia becomes less and less significant so that the PSS approximation can be considered as reference (we remind that Neumann's solution is computed with a semi-infinite domain so that no steady state can be reached in this case).
As a conclusion, the numerical simulations appear quite satisfactory since they stick to Neumann's solution in the initial time steps and to the PSS solution in later stages of the process.

Applications to polar-symmetric problems
Below, some XFEM applications to problems with cylindrical and spherical symmetry are presented.The one-dimensional XFEM model used in the previous Sections can be adapted to the present cases by considering the specific symmetry and suitably implanting the volume integration in Eqs (37).The one-phase melting of a cylinder and a sphere of iced water are considered in Sections 7.1 and 7.3, respectively.Since no exact solutions exist for these cases, a sufficiently large value of Stefan number is used, so that the PSS approximations, which are developed in the relevant Sections, can be considered for comparing the numerical results.In Section 7.2, the problem of technical interest of the freezing water into a thermally-insulated pipe is studied.

Melting of iced-water cylinder
We study the melting of an iced water cylinder, whose cross-section is shown in Figure 21.The cylinder is initially at the melting temperature and an essential boundary condition   1 t r a n

T R T T
  is imposed to the external surface so that the initially solid cylinder progressively melts.
The problem has cylindrical symmetric and can be described by the spatial coordinate r .The energy balance equation in cylindrical coordinates and for isotropic heat conduction is The phase-interface temperature condition is The boundary conditions are We develop the PSS approximation to compare the numerical results.Under the PSS hypotheses, Eq. ( 65) becomes By integrating and applying the boundary conditions of Eqs (68), the approximated temperature field is obtained as By substituting in Eq. (66) and integrating it with the initial condition in Eqs (69), an implicit relation between the interface position and the time is obtained, For the numerical analysis we consider the data in Table 1, from which a relatively high value of Stefan number is achieved, so that the PSS approximation represents an adequate reference for comparing the numerical results.
Figure 22 shows the interface position calculated through the PSS approximation and the XFEM simulations with two different spatial discretisations, with 10 and 40 equally-sized finite elements.The increase in the spatial discretisation improves the quality of the numerical solution especially in the last stages of the melting process, when the phase-interface is close to the cylinder axis.This is mainly due to the evaluation of the integrals in Eqs (37).In particular, the accuracy of a finite element is higher if the areas of the cylindrical surfaces constructed with the radii at the two nodes of the considered element are closer.The temperature field at some time instants, 2400 s, 4800 s ,…, 12000 s t  , is represented in Figure 23.Temperature fields evaluated through PSS approximation and XFEM simulation differ, if they are evaluated at the same time instant, because of the error in the phase-interface position.However, the temperature fields superimpose very well, if they are 'artificially' compared for the same phase-interface position (at different time instants), as shown in Figure 24, attesting the correctness of the simulated temperature profile.

Water freezing in thermally-insulated pipe
The solidification of water in a thin metal pipe is considered.The water is assumed to be at rest and the pipe radium is taken to be 0.05 m.We assume the pipe to be surrounded by a 0.01-m-thick thermally-insulating layer as depicted in Figure 25.The initial temperature of the water is taken to be equal to the freezing temperature and the ambient temperature equal to 10 °C  (this results in an essential boundary condition at 0.06 m r  ).Further data of the problem are reported in Table 1, where the subscript 'ins' denotes the thermally-insulating material.The presence of the metal pipe is neglected due to its small thickness and its high thermal conductivity.The XFEM model is constructed by subdividing the ice/water domain into 5 finite elements and by using an additional element for the thermally-insulating material.The time discretisation used is  1 10 s, 2 10 s, 3 10 s t     .Within the considered domain, there are two gradient discontinuities in the temperature field; the first one is due to the moving phase interface and the second one, to the different thermal conductivity in the ice and insulating material.The evolution of the interface position within the pipe is plotted in Figure 27.For this example no experimental or analytical comparison is available.However, it is reasonable to consider the present simulations as fairly close to the physical case, as the results presented in the previous Sections show that the numerical model behaves very well in a variety of planar-symmetric examples (Section 6), and it has been validated for cylindrical symmetry (Section 7.1).

Melting of iced-water sphere
Moulds for ice spheres are advertised for serious on-the-rock drinkers by asserting that the minimum surface-to-volume ratio means that the ice melts at a slower pace, preventing the drinks from getting watery too quickly.This seems questionable and, indeed, the problem is not well posed, as the amount of melted ice depends mainly on the heat transferred from the container to the environment.Here, we present a simpler problem, where a sphere of water ice melts in liquid water.We assume the heat transferred by convection to be negligible.
The present problem with spherical symmetry, schematised in Figure 28, is analogous to the cylindrical one presented in Section 7.1.The sphere is initially at the melting temperature and the essential boundary condition   1 t r a n T R T T   is imposed, so that the initially solid sphere starts melting.Due to the spherical symmetry of the problem, it can be studied by considering only the coordinate r (we use the same symbol r for the radial coordinate in both cylindrical and spherical coordinate systems, the context ensuring the disambiguation).The energy balance in spherical coordinates and for isotropic heat conduction is whereas Stefan law, phase-interface temperature condition, boundary and initial conditions read as in Eqs from (66) to (69).Since no analytical solution is available.We develop then the PSS approximation.Eq. (74) becomes For the numerical analysis, the same data used in the cylindrical example in Section 7.1 are adopted (see Table 1), so that the same value of Stefan number is obtained, high enough for considering the PSS approximation as a suitable reference for comparing the numerical results.Figure 29 shows the phase-interface position calculated through the PSS approximation and the XFEM simulations with 10-and 40-element discretisations.As for the cylindrical problem discussed in Section 7.1, the approximation with the finer spatial discretisation improves especially in the vicinity of the sphere centre.In this case, such behaviour is even more pronounced, as it can be evinced by comparing Figure 29 and Figure 22.In fact, the accuracy depends on the ratios between the areas of the surfaces (cylindrical or spherical) at the two nodes of one element; the cylindrical case has then a better behaviour than the spherical one, as the areas are proportional to the radii, instead of the squares of the radii.The temperature field at 2400 s, 4800 s, 7200 s, 9600 s t  is shown in Figure 30.Temperature fields evaluated through PSS approximation and XFEM simulation notably differ, if they are evaluated at the same time instant, because of the error in the phase-interface position.Still, the temperature fields compare very well, if they are evaluated for the same phase-interface position (but at different time instants), as shown in Figure 31.

Remarks
The Stefan problem for melting/solidification has been approached through two numerical techniques that enhance the FEM by allowing a temperature gradient discontinuity across the phase interface without remeshing.In the XFEM, the discontinuity is introduced through additional nodal shape functions, whereas in the embedded-discontinuity approach it is controlled at elemental level.In both cases the position and evolution of the phase interface can be described through the level set method.The presented analyses are restricted to one-dimensional cases; therefore, only a speculative discussion can be offered for higher-dimensional problems, which are in fact object of intensive research at the present time.
As to the XFEM, the adopted approach is the one proposed by Chessa et al. [16], for which, here, different numerical techniques for the constraint enforcement have been used and compared; furthermore, a detailed numerical study of the error and a comparative review have been performed.Though the penalty approach can be considered for practical applications, Lagrange multiplier method enforces the constraint exactly and does not require the calibration of an unphysical penalty coefficient.Lagrange multiplier method can be directly applied to the discretised constraint; however, as it is inferred from the quoted references, in higher dimensional cases attention should be paid in the choice of Lagrange multiplier space in order to achieve stability.The interpolation to be adopted for the blending elements has been also discussed; it is concluded that if no enrichment is used in these elements the accuracy is improved without any particular disadvantage; this should extend also to higher dimensional cases.
As to the proposed embedded-discontinuity method, it has been shown to have a comparable accuracy and simplicity of implementation with respect to the XFEM.The generalisation of this approach to higher-dimensional problems is possible by constructing suitable shape functions; however, this is less straightforward than in the case of the XFEM, where the enrichment functions can be constructed elegantly and automatically by using Eq.(32).

Figure 1 .
Figure 1.Schematic representation of the domain.


The domain  is partitioned into a set of finite elements   e

Figure 2 .Figure 3 .
Figure 2. Support of a standard FE shape function and its isocontours in a 2D bilinear case.

1 
the case of constant c .They perform an error study by testing different values of the parameter the highest accuracy is obtained for  , i.e. in the case of backward Euler algorithm.

Figure 4 .
Figure 4. Shape functions within one element in local coordinate x: (a) standard, (b) enriched in XFEM, and (c) enriched in embedded-discontinuity FEM.

Figure 5 .
Figure 5.Standard FE shape functions and embedded-discontinuity function.

Figure 8 .
Figure 8. Convergence of the interface position.

Figure 9 .
Figure 9. Relative error in the interface position.

Figure 10 .
Figure 10.Error in temperature field versus element size.

Figure 12 .
Figure 12.Temperature field in the vicinity of the interface.

Figure 13 .Figure 14 .
Figure 13.Error versus element size with different kind of blending elements.
analysis, we discretise the spatial domain into 10 equally-sized elements and time with

Figure 15 .
Figure 15.Temperature field in the element containing the discontinuity.

Figure 17 .
Figure 17.Error in temperature field versus element size for embedded-discontinuity FEM.

.
The boundary conditions are   0 t r a n 0, T t T T   and  1 t r a n , T t T T   

Figure 20 .
Figure 20.Temperature field at different time instants.
and the position of the phase interface   * r t is governed by Stefan law,

Figure 21 .
Figure 21.Cross-section of the melting ice cylinder.

Figure 22 .
Figure 22.Interface position in the cylindrical melting example.

Figure 27 .
Figure 27.Interface position within the pipe.
integrating and considering the boundary conditions, the temperature field obtained.Substituting in Eq. (66) and integrating it with the initial condition of Eqs (69) leads to the implicit relation for the time evolution of the phase-interface,

Figure 30 .
Figure 30.Temperature field within the spherical domain.

Figure 31 .
Figure 31.Temperature field for a given phase-interface position in spherical symmetry.

Table 1 .
is initially at the freezing/melting temperature, with Stefan number is computed with the temperature difference between phase-transition and coolant thermostat,

Table 1 .
Data used in the numerical examples.