## 6 Numerical study on the error and applications to planar-symmetric problems
### 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-17], because of its robustness and simplicity of implementation. As discussed in  
Figure 10: Error in temperature field versus element size.  
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  
\[\varepsilon_{\text{\tiny{constr}}}\left(t\right)=\frac{\left|T_{\text{\tiny{ tran}}}-T\left(x_{\text{\tiny{c}}}\left(t\right)\right)\right|}{\left|T_{\text{\tiny{ tran}}}-T_{0}\right|}\,. \tag{58}\]  
The conditioning of the linear system in Eq. (42) is measured by the condition number with respect to inversion, \(\alpha_{\text{\tiny{cond}}}\), defined as the ratio between the maximum and the minimum singular values of the system matrix \(\vec{\mathbf{K}}+b\mathbf{K}_{\text{\tiny{pen}}}\) ; the higher is \(\alpha_{\text{\tiny{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  
\[\beta=b\left(\frac{k_{\text{\tiny{sol}}}+k_{\text{\tiny{fig}}}}{2}\frac{S}{ \ell_{\text{\tiny{c}}}}\right)^{-1}\,, \tag{59}\]  
Figure 11: Penalty method: relative error **(a)** and condition number **(b)**.  
with \(S\) and \(\ell_{e}\) the average cross-sectional area and length of the one-dimensional elements, respectively. Since \(\mathbf{K}_{\text{pen}}\) is dimensionless, the term in parentheses in Eq. (59) contains the physical dimension and the average magnitude of the coefficients of the stiffness matrix \(\mathbf{K}\), which, in turn, contributes to \(\mathbf{\tilde{K}}\). In this way, \(\beta\) 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 \(\beta\) 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 \(\beta\) 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 \(t\) = 6\(\times\)10\({}^{5}\) s as obtained with penalty ( \(\beta\) = 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 \(\beta\) = 18, with an error in the constrain enforcement exceeding 5%.