A XFEM LAGRANGE MULTIPLIER TECHNIQUE FOR STEFAN PROBLEMS

The two dimensional phase change problem was solved using the extended finite element method with a Lagrange formulation to apply the interface boundary condition. The Lagrange multiplier space is identical to the solution space and does not require stabilization. The solid-liquid interface velocity is determined by the jump in heat flux across the i nterface. Two methods to calculate the jump are used and compared. The first is based on an averaged temperature gradient near the interface. The second uses the Lagrange multiplier values to evaluate the jump. The Lagrange multiplier based approach was shown to be more robust and precise.


INTRODUCTION
The finite element method Reddy (2006) has been extensively studied and successfully used in a wide variety of scenarios involving continuous media but particular situations are still problematic. The finite element method uses polynomial interpolations within individual elements to approximate the solution. Consequently, it can only be applied to problems with discontinuities by splitting the domain into submeshes. This makes the finite element method ill suited to solve problems involving discontinuities that are part of the solution or move in time. The Stefan problem Nedjar (2002); Beckermann et al. (1999); Helenbrook (2013); Özişik (1993) for the isothermal solidification or melting of a material is one such situation.
The extended finite element method Belytschko et al. (2001); Dolbow et al. (2000); Belytschko et al. (2009) is based on the partition of unity method Babuska and Melenk (1997); Dolbow et al. (2000); Melenk and Babuska (1996). Using carefully selected functions ψ(x, t), the technique adds additional degrees of freedom that will "enrich" the interpolation and allow the solution to adopt a discontinuous behavior. The particular type of behavior is determined by the enrichment function ψ(x, t), known a priori. Only those nodes who's support is cut by the interface and have a modified behavior must be enriched (see figure 1), making the additional computational costs local to the interface. The interface geometry is stored and transported in a computationally efficient manner, most commonly using the level set method Osher and Sethian (1988); Osher and Fedkiw (2001).
An important challenge in the extended finite element method is the imposition of Dirichlet boundary conditions on the interface. The absence of nodes on the interface means that we cannot apply specified values directly and an additional constraint must be added in the finite element formulation to apply the appropriate boundary condition on the interface. The two main numerical techniques used are the penalty method Chessa et al. (2002); Zabaras et al. (2006); Bernauer and Herzog (2011) and Lagrange multiplier method Moes et al. (2006); Ji and Dolbow (2004); Babuska (1973). The penalty method requires the definition of a free numerical parameter to be determined by trial and error. The Lagrange multiplier requires no such numerical parameter but is more computational expensive and may present oscillations in the solution near the interface if an improper interpolation space is used for the multiplier Moes et al. (2006); Ji and Dolbow (2004).
A recent effort has been done by Gerstenberger and Wall Gerstenberger and Wall (2010) to eliminate this obstacle and facilitate the use of Lagrange multipliers. They have developed a Lagrange multiplier formulation for the solution of problems involving voids in the geometry. The main advantage of their approach is the use of identical interpolation spaces for the solution and Lagrange fields. Although no mathematical proof of its stability was given, numerical applications with the stationary diffusion equation Gerstenberger (2010) and Navier-Stokes equations in fluid-structure interactions Gerstenberger and Wall (2010) were shown to be stable.
In the work presented here, this novel Lagrange multiplier technique is applied to the classical two dimensional Stefan phase-change problem. The physical nature of the problem being quite different to the FSI problem Gerstenberger and Wall (2008), a different type of enrichment scheme is required. Specifically, a weakly discontinuous temperature field and strongly discontinuous Lagrange multiplier field at the liquid/solid interface. The jump in the Lagrange multiplier value at the interface is then used to determine the interface velocity and compared with a temperature gradient based heat flux calculation.
The paper is divided as follows. The governing equations for the Stefan problem are described in section 2 . The finite element formulation, level set problem and details concerning the interface movement and extended finite element method are described in section 3. Benchmark examples are then solved in section 4 to validate the new Lagrange multiplier approach and compare its performance with the penalization technique commonly found in the literature. Finally, the paper ends with some concluding remarks in section 5.

Problem Formulation
Consider a domain Ω with an initial temperature T (x, t0) and interface Γ separating solid (Ωs) and liquid (Ω l ) phases with different thermal properties. We suppose that the density is identical in both phases and that the material has an isothermal phase change at some melting temperature Tm. Applying the conservation of energy in Ω results in equation (1a), where (cp)i, i = l, s is the specific heat, ki, i = l, s the thermal conductivity and ρ the density. Additionally, the melting temperature must be applied on the solid-liquid interface (1b). Dirichlet and Neumann type boundaries away from the interface are applied on ∂Ω = ΓN ∪ ΓD as usual (1c,1d).
Conservation of energy at the interface requires that the jump in heat flux normal to the interface (caused by the imposition of the melting temperature) be related to the rate of solidification or melting of the material as described in equation (2), where L is the latent heat and vΓ the normal interface velocity Özişik (1993).
The normal vector ns points from the liquid to solid phase, meaning that the interface velocity is positive for melting and negative for solidification.
Tracking the moving interface is done using the level set method Osher and Fedkiw (2001); Osher and Sethian (1988). The principle behind this method is to introduce a new variable φ(x, t) defined as the signed distance function to the interface (3). The interface is then easily identified as the set of points where φ(x, t) = 0.

Enriched Interpolation Scheme
To account for the jump in heat flux at the interface, the temperature field must be continuous with a jump in the gradient. This behavior is captured by using the approximation (4) for the temperature field Chessa et al. (2002) where Ni are the standard interpolation functions, Ti and T * j the standard and enriched degrees of freedom, respectively and ψj(x, t) the enrichment function, based on the absolute value of the level set field.
When using (4) special attention must be given to elements containing enriched nodes that are not cut by the interface, called blending elements. A modified interpolation scheme must be used in these elements to maintain an optimal convergence rate, as described in Fries (2008); Shibanuma and Utsunomiya (2009). A more compact way to write (4) is by use the more standard matrix form (5).
The Lagrange multiplier q used to apply the fusion temperature, as described in (13), must reproduce the behavior of the heat flux, which will have a jump at the interface. This behavior is captured by using the approximation (6). In this case, the enrichment function is constructed with a modified Heaviside function.
Following the expression (5), the Lagrange multiplier may be written as (7) where the symbol is used to distinguish between the temperature and Lagrange field interpolation functions.

Finite Element Formulation
The weak form of the energy conservation equation (1a) is (8), where δT is the test function evaluated at the current time step: δT = δT (t n+1 ) Fries and Zilian (2009).
Using a backward Euler scheme for the time derivative of T in (8) gives Fries and Zilian (2009): Substituting the approximation for the temperature field into (9) leads to the finite element system of equations (10).
In elements which are intersected by the interface, an additional constraint must be applied to the formulation to take into account the interface boundary condition (1b). In this work, two methods are used. The first is the penalization method Chessa et al. (2002); Bernauer and Herzog (2011);Zabaras et al. (2006), which applies the melting temperature on the interface by multiplying (1b) by a very large penalization parameter β (eq. 11) and including it in the weak form (8) for intersected elements only.
The complete system of equations for intersected element then becomes: This method is simple to implement and adds very little computational effort. However, the choice of β can be an important factor in the solution's precision. If β is too small, the constraint will not be properly taken into account. If β is too large, oscillations can appear along the interface Ji and Dolbow (2004). Moreover, the optimal value is problem dependent and must be found by trial and error.
The second method used in this work to impose the proper boundary condition on the interface is the Lagrange multiplier Ji and Dolbow (2004). This method adds a secondary variable to the formulation. Physically, this secondary variable corresponds to the heat flux generated on the interface from the additional constraint on the problem. Normally, this secondary flux variable is defined purely on the interface and requires the primary (temperature) and secondary (heat flux) variables to respect the inf-sup condition Babuska (1969); Brezzi (1974). Otherwise, oscillations may appear in the solution near the interface Béchet et al. (2009); Ji and Dolbow (2004); Moes et al. (2006).
To overcome this difficulty, a new adaptation of this method was developed in Gerstenberger and Wall (2010); Baiges et al. (2012), based on the Lagrange multiplier technique found in Zilian and Fries (2009). Here, the Lagrange multiplier is defined as a vectorial flux and interpolated on the same mesh as the temperature field. The projection of this secondary variable on the interface is then used as a scalar Lagrange multiplier to impose the melting temperature. Note that since the secondary variable is now a vector field, two additional unknowns have been added. To obtain a properly defined problem, the secondary variable is weakly coupled with the flux calculated from the temperature gradient in the domain and a complete system of equations (13) is obtained Gerstenberger (2010).
After applying the backward Euler scheme Fries and Zilian (2009) to the time derivative and replacing T and q with their interpolation schemes (5) and (7), we obtain the following system of equations for intersected elements: In elements which are not cut by the interface, the Lagrange multiplier is weakly coupled with the temperature gradient but no constraint is present and the system reduces to: As described in Gerstenberger and Wall (2010), the interpolation used for the Lagrange multiplier can be C -1 continuous (continuous temperature, discontinuous heat flux) at inter-element boundaries, allowing the condensation of equation (14) on the element level. The resulting contribution of q is added to the global matrix for elements intersected by the interface. Only the temperature field is solved, reducing the size of the global system of equations. The secondary flux variables may then be calculated from the temperature values Gerstenberger and Wall (2010).

Level Set Formulation
Once an initial value φ(x, t0) is defined using (3), the interface movement is governed by its transport equation (16a) where v is the convection velocity and F (16b) is the interface speed in the normal direction. The calculation of F is explained below.
Equation (16a) is solved explicitly (forward Euler scheme) using a linear interpolation. The weak form is given by: In most applications, the normal component F is only known on Γ. In order to solve (17) on Ω, a valid value for F must first be constructed using the following problem Chessa et al. (2002): This approach guarantees that the φ field velocity is everywhere normal to the interface and is coherent with the interface's physically determined velocity. For more details concerning the construction of F, see Osher and Fedkiw (2003); Chessa et al. (2002). In this paper, two different methods are used to evalute vΓ and are described below.
Equation (16a) is first order hyperbolic and must be stabilized to minimize the presence of oscillations in the solution Chessa et al. (2002); Bernauer and Herzog (2011). The GLS method is used here Hughes et al. (1989). The level set method offers several advantages. It is easily extensible to three dimensions and stores the interface location as a scalar variable. Furthermore, the level set field can be defined in a small region surrounding the interface and the level set formulation solved locally, reducing the impact on the total simulation computation time. It is also robust enough to handle interface merging and breaking naturally Osher and Fedkiw (2001).
The main disadvantage of the the level set method is its tendency to deviate from a signed distance function over time Osher and Fedkiw (2001). This error accumulates with additional time steps and degrades the quality of the solution, particularly the level set gradient near the interface. This distortion can be a source of error in the numerical solution of the level set formulation and the physical problem on which it is based. Therefore, it is necessary to reinitialize φ(x, t) regularly to maintain an acceptable solution ( ∇φ ≈ 1). Another limitation to the level set method is the use of an explicit time scheme, which limits the size of the time step. The explicit time step is required in order to determine the nodes to enrich. In other words, the interface position must be determined before equations (12) or (14) can be solved.

Interface velocity calculation
The proper evaluation of fluxes on either side of the interface in equation (2) is crucial in obtaining a precise and robust model. The simplest way to evaluate the jump in heat flux across the discontinuity is to evaluate the gradient at appropriately chosen points in the solid and liquid phases at some small distance away from the interface Chessa et al. (2002); Zabaras et al. (2006). However, this approach may be the least accurate option Ji and Dolbow (2004). A more involved but robust technique is to evaluate the temperature at multiple points at specific distances from the interface to obtained an averaged value Bernauer and Herzog (2011). This approach is used in this work and given in (19), where δx is some fraction of the average element size (figure 2). vΓ = 1 ρL When applying the Lagrange multiplier technique another option becomes available. As previously mentioned, the Lagrange multiplier corresponds to the heat flux within the element. We can evaluate the jump in  (20).
This is done by evaluating the Lagrange field on the interface, approaching from either side. A similar strategy was used by Merle and Dolbow (2002), using an iterative procedure based on the LATIN method to impose the interface temperature. To our knowledge, no Lagrange multiplier based algorithm for the Stefan problem has been developed using this strategy for the interface velocity. The final algorithm can be described as follows. Assuming a given time t n , temperature solution T n and level set solution φ n , the strategy to solve for T n+1 consists in the following steps: 1. Compute the interface velocity v n Γ using (19) or (20) 2. Construct F on the level set domain by solving problem (18) 3. Solve for φ n+1 using (17) 4. Solve for T n+1 using (12) or (14) 5. Set t n = t n+1 and go to step 1

Numerical Integration
The introduction of discontinuous functions inside elements greatly reduces the precision of standard Gaussian quadrature and may lead to rank deficient matrices âȂłChessa et al. (2002). An accurate but geometrically complex solution is to subdivide elements involving discontinuities into continuous subelements Moes et al. (1999); Chessa et al. (2002); Gerstenberger and Wall (2010). Each element is subdivided into a number of subelements (lines, triangles or tetrahedrons), as shown in figure 3, to properly fit the contour of the interface (point, line or surface) and element boundaries. The integral over the entire element Ie is then the sum of the integration of each subelement Is using standard Hammer quadrature. It is important to note that subelements carry no degrees of freedom or interpolation functions. They are only required as a geometrical tool to construct the element integrals. In transient problems the location of the quadrature points must change as the interface moves in time, requiring that every cut element be subdivided at each time step. However, the subdivision is applied only to a small number of elements, reducing the overall increase in computational effort required.
In transient problems, the interpolation functions at time steps n and n + 1 are based on different positions of the interface and are discontinuous at different places in the element. The integration scheme for the mass matrix (equation (10c)) must take both intersections into account when generating the integration subelements to obtain optimal convergence Fries and Zilian (2009). This can be difficult and can significantly increases the number of subelements required to fit the geometry. However, previous authors have successfully used integration schemes considering the current interface position only Chessa et al. (2002); Chessa and Belytschko (2003) and this strategy is used in this work. As suggested in Fries and Zilian (2009), the test functions are evaluated using the current time step's level set values.

RESULTS
To validate the new Lagrange formulation and compare its performance with the penalty method, two benchmark problems were solved with both approaches. To evaluate the precision of the imposed Dirichlet boundary condition on the interface, the penalty (12) and Lagrange formulations (14) were compared using the same gradient based interface velocity algorithm (19), referred to as case I. A third solution to the problem was also obtained by using the Lagrange multiplier jump value to evaluated the interface velocity (20), referred to as case II .
To evaluate the impact of the Lagrange field's polynomial degree on the solution, two implementations of the Lagrange formulation were used. In the first, the Lagrange field is linear and continuous at interelement boundaries (C 0 continuous). In the second, the Lagrange field is constant per element and condensed on the element level (C -1 continuous) Gerstenberger and Wall (2010). Both linear and constant interpolation schemes were tested with both gradient (19) and Lagrange (20) based velocity values, for a total of 4 distinct solving algorithms based on the Lagrange multiplier technique and one using penalization. The temperature field interpolation is linear in all cases.

2 Phase 1D problem
The first benchmark problem is the one dimensional two phase analytical solution of the Stefan problem in a semi-infinite domain (x > 0), taken from Merle and Dolbow (2002). The thermal properties are phase dependent and given in table 1. The domain is initially liquid, as shown in figure 4, with the solid-liquid interface 2 mm away from the left domain boundary. The initial temperature is 277 K. The top and bottom edges are insulated. At t = 0, the temperature on the left edge is lowered to 263 K and the right edge is maintained at 277 K. Temperature evaluation points for the gradient based velocity calculations are taken at a maximum distance of 10% of the mean element size and β = 10 8 where applicable. The interface position, as a function of time for both interface boundary condition techniques, is shown in figure 5. Three different mesh sizes were used to verify the convergence of the various techniques, and are shown in figure 6. The error norms are calculated using (21a) where x a Γ is the analytical interface position and λ = 0.3073 Merle and Dolbow (2002). The time step was chosen as ∆t = ∆x 2 ρcp k s Merle and Dolbow (2002).
As shown in figure 5, the numerical solutions follows the analytical solution up to approximately 3000 seconds. Beyond this point, the effects of the finite computation domain become apparent. The same behavior was observed by Merle and Dolbow (2002). Consequently, the  To obtain a precise value, the flux jump must be calculated directly on the interface. The use of the Lagrange multiplier gives us an easy and precise way of doing so. The use of different interpolation schemes has a much smaller impact on performance. A small increase in precision is obtained by using a linear interpolation for the Lagrange multiplier field, no matter the interface velocity used. In fact, the use of a constant per element Lagrange field leads to a slightly less precise solution than the penalization technique when the gradient based interface velocity is used. This difference may be caused by the use of bi-linear quadrilateral elements for the temperature field. When using bi-linear interpolation functions, the gradient varies linearly inside the element. Consequently, the constant per element Lagrange field cannot reproduced the behavior of the temperature gradient within the element exactly. An identical performance is observed between the linear Lagrange field and penalization algorithms using the gradient based interface velocity, indicating that both approaches enforce the melting temperature equivalently.
From a more practical point of view, the gradient based calculation has certain drawbacks. First of all, it requires the determination of 6 evaluation points normal to the interface, which will vary with changing interface geometries and must be determined at every time step. Secondly, some of theses points may be outside of the calculation domain depending on the interface position (see figure 2) and must be treated by some approximation, depending on the situation. These extra tasks increase the computational effort and code complexity. By contrast, the Lagrange based velocity is evaluated directly on the interface and avoids these extra steps and approximations. Ultimately, the extra computational effort required for the gradient based interface velocity is traded for the effort required to solve the more involved condensed Lagrange formulation. Depending on the exact implementation, either approach may be slightly faster. In our case, the condensed Lagrange formulation was slightly quicker. The uncondensed Lagrange formulation lead to a much larger tangent matrix and was the slowest due to the longer global system resolution time.
Time steps significantly smaller than ∆x 2 ρcp k s would lead to erroneous interface temperatures at early times with the Lagrange formulation as shown in figure 7. The effect of time step size on the solution of transient heat transfer problems using the finite element method has been previously investigated in Thomas and Zhou (1997). However, it is important to note that the penalty formulation was unaffected at the tested time step of ∆x 2 100 ρcp k s (figure 7), suggesting that the Lagrange multiplier is more sensitive to this phenomenon.

2 Phase 2D problem
The second benchmark problem is the two phase analytical solution of the Stefan problem in two dimensions, taken from Aysoufi and Keith (2003). The thermal properties are constant and given in table 2. The domain is initially liquid, as shown in figure 8. The initial temperature is 273.3 K. The top and right edges are insulated. At t = 0, the temperature on the left and lower edges are lowered to 272 K. The analytical solution of this problem was first developed in Rathjen and Jiji (1971). The nondimensionalized interface position y (x ) is determined using equation (22c) where C = 0.159, m = 5.02 and λ = 0.70766 for the given material properties and boundary conditions, α is the thermal diffusivity and x the non-dimensionalized x axis.  Temperature evaluation points for the gradient based velocity calculations are taken at a maximum distance of 35% of the mean element size and β = 10 9 where applicable. The level set formulation used in this work does not have a reinitialization procedure. To validate the present Lagrange formulation and compare the various algorithms, the problem was simulated using ∆t =5 × 10 −5 s. The final interface position at t = 0.025s is shown for the different algorithms and the analytical solution in figure 11. The error norms for the final interface position are given in table 3. The error norms were calculated using equation (23) where (x num n , y num n ) is the position on the interface of the numerical solution, (x min , y min ) the position on the analytical interface (22a) closest to (x num n , y num n ) and n the total number of points taken on the interface.
As shown in figure 11, the numerical solution is in agreement with the analytical solution for all algorithms, indicating that both approaches  Fig. 9 Interface position error norm for 2D problem enforce the melting temperature appropriately in 2D as well. No improvement in the solution is obtain by using the linear or constant Lagrange formulation compared to the penalty formulation when using the gradient based interface velocity. Using the linear Lagrange field interface velocity significantly improved the solution. The use of the constant Lagrange field interface velocity also increased the accuracy but to a smaller extent. This is probably due to the inability of the constant Lagrange field to reproduce the linear gradient of the bi-linear quadrangle interpolation used for the temperature. As a result, a linear Lagrange formulation, using a Lagrange based interface velocity, would lead to an optimal algorithm in terms of accuracy. To reduce the resolution time, a C -1 continuous linear interpolation could be used and condensed on the element level Gerstenberger and Wall (2010). It is well known in the literature that the level set field must be regularly reinitialized to maintain accurate results Chessa et al. (2002); Osher and Fedkiw (2003). The one dimensional problem does not require this step because the interface velocity is constant on the entire domain. The two dimensional problem involves varying interface velocities in space and distortions in the level set field may appear. The absence of a reinitialization step in this work did not have a significant impact on the stability of the simulation.
A detailed convergence study for different time steps and mesh sizes would be influenced by the accumulated error in the level set field, making comparisons between algorithms difficult. However, a convergence study was done using the linear Lagrange field formulation for completeness. The results are shown in figure 9 for the element size and figure 10 for the time step. As we can see, the Lagrange formulation converges as the mesh is refined, although the convergence rate appears to be lower than in the one dimensional case. Decreasing the time step size has the reverse affect, increasing the error in the solution. These sub-optimal convergence results are, at least in part, caused by the absence of the level set reinitialization. As the simulation advances, errors in the level set field accumulate, leading to errors in the overall solution. Decreasing the time step size increases the amount of time steps needed to reach the final time, quickening the degrading of the level set solution and increasing the error. The one dimensional problem is immune to the effect because the linear interface creates a uniform level set velocity field and no distortions can appear.
Another practical aspect that was observed in the 2D case was the importance of numerical parameter selection. For the penalty formulation, two user-based numerical parameters are required: the penalty term β and the maximum element distance δx to evaluate the interface velocity. An inappropriate selection of either of these parameters could lead to a significant reduction in precision of the interface position, with the dis-  Fig. 11 Final non-dimensionalized interface position for 2D problem tance δx having a much more significant impact than β. As an example, the benchmark problem using the penalty formulation was solved using δx = 20%. The interface position at t = 6.25 × 10 −3 s, at which point the program fails, is shown in figure 12.

CONCLUSION
In this work, the Lagrange multiplier scheme developed by Gerstenberger and Wall (2010) was applied to the classical Stefan problem in one and two dimensions using a weakly discontinuous temperature field and strongly discontinuous flux-based Lagrange field. We have shown that it allows the application of Dirichlet boundary conditions with comparable precision to the penalty technique Chessa et al. (2002) using constant and linear interpolation schemes for the Lagrange field. Furthermore, a gain in accuracy was obtained by calculating the interface velocity using the La-  Fig. 12 Interface position at t =6.25 × 10 −3 using δx = 20% grange field instead of the more widely used temperature gradient and reduced the discrepancy in computational cost between the two techniques. The Lagrange based formulation requires no user defined numerical parameters, improving the model's robustness. Further work will be realized to account for convection in the liquid phase and mass flux at the interface in problems involving different densities in the solid and liquid phases. the Fonds de recherche du QuÃl'bec -Nature et Technologies by the intermediary of the Aluminium Research Centre âȂŞ REGAL. Special thanks to Patrice Goulet for his invaluable assistance with the development of the level set algorithm, and advice on the art of code writing.