In the second appr~ach~,~ a deforming spatial mesh is used and a three-stage iterative cycle isfollowed. First, a shape of the free boundary is assumed. Secondly, the boundary valueproblem on the resulting domain is solved after disregarding one of the boundary conditionsat the free boundary. Thirdly, the shape of the free boundary is updated using the previouslyneglected boundary condition. This iterative cycle is repeated until some desired convergencecriterion is satisfied. This approach suffers from a number of disadvantages. The iterativecycle often does not converge, and even when it does, being a fixed point type iteration, itconverges slowly. Moreover, a new finite element problem must be solved at each iteration.(3) The third, more recent approa~h~-~ also involves a deforming mesh but eliminates thesuccessive iteration between the free boundary position and the field variables by introducingthe position of the nodes at the free boundary directly as degrees of freedom. The non-linearequations are then solved using a Newton-Raphson (or a quasi-Newtori9) iterativeprocedure which results in the simultaneous calculation of the position of the free boundaryand the field variables. The advantage of this method is its second-order (or superlinear forthe quasi-Newton method) rate of convergence. A disadvantage is that a complete account ofthe variations with respect to the free boundary degrees of freedom must be incorporated intothe Jacobian of the system of equations. These variations, which involve integrals over a largepart of the domain, have a non-local character. This means that the method does not fit intostandard finite element codes, where the coefficient in the equations of an unknownbelonging to a nodal point is determined completely by contributions over the neighbouringelements. This makes the implementation of the method relatively difficult. Anotherdisadvantage is that a new finite element problem must be solved at each iteration.The aim of this paper is to derive a numerical method, the total linearization method (TLM),which is much easier to implement than the Newton-Raphson algorithm while retaining itssuperior convergence properties. It is similar to References 10 and 11, but in this paper it will beshown that the influence of the unknown position of the free boundary can be reduced completelyto boundary integrals. This has great advantages for software implementation.The sequel is restricted to the die-swell problem, which serves as a model case for the TLM.However, it is possible to apply the TLM to other free boundary problems as well.In Section 2 the die-swell problem is described and its mathematical formulation is given.In Section 3 the weak formulation of the problem is derived. This weak formulation is linearized inSection 4. The numerical method consists of a discretization of this linearized weak formulation. InSection 5 the numerical results are shown and compared with those available from the literature.