# Phase-field and lip-field approaches for fracture with extreme mesh deformation (X-Mesh): a one-dimensional study

Nicolas Moës, Benoît Lé, Nicolas Chevaugeon, Jean-François Remacle

► To cite this version:

Nicolas Moës, Benoît Lé, Nicolas Chevaugeon, Jean-François Remacle. Phase-field and lip-field approaches for fracture with extreme mesh deformation (X-Mesh): a one-dimensional study. Comptes Rendus. Mécanique, 2025, 10.5802/crmeca.318. hal-05262622

HAL Id: hal-05262622

https://hal.science/hal-05262622v1

Submitted on 16 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Comptes Rendus
## Mécanique

Nicolas Moës, Benoît Lé, Nicolas Chevaugeon and Jean-François Remacle

**Phase-field and lip-field approaches for fracture with extreme mesh deformation (X-Mesh): a one-dimensional study**

Volume 353 (2025), p. 1027-1051

Online since: 16 September 2025

https://doi.org/10.5802/crmeca.318

<figure><img src="image_2.png" /></figure>

This article is licensed under the

CREATIVE COMMONS ATTRIBUTION 4.0 INTERNATIONAL LICENSE.

http://creativecommons.org/licenses/by/4.0/
Research article / *Article de recherche*

# Phase-field and lip-field approaches for fracture with extreme mesh deformation (X-Mesh): a one-dimensional study

## *Approches phase-field et lip-field pour la rupture avec déformation extrême du maillage (X-Mesh) : une étude unidimensionnelle*

Nicolas Moës<sup>*</sup>,<sup>a,∗</sup>, Benoît Lé<sup>b</sup>, Nicolas Chevaugeon<sup>c</sup> and Jean-François Remacle<sup>d</sup>

**Abstract.** We consider a one-dimensional fracture problem modelled using either the phase-field or lip-field approach. In both cases, we optimise the incremental potential with respect to the displacement and damage fields and the nodal coordinates of the mesh. This is thus a variational mesh study. We observe that, as the damage reaches its maximum value, the optimisation drives the most damaged element to zero size as the damage reaches its maximum value. This peculiar element provides a precise displacement jump representation as the bar breaks. The overall solution is also shown to be much more accurate than the fixed mesh solution. This work forms part of an exploration into the capabilities of extreme meshes in computational mechanics (X-Mesh).

**Résumé.** Nous considérons un problème unidimensionnel de rupture modélisé à l'aide de l'approche *phase-field* ou *lip-field*. Dans les deux cas, nous optimisons le potentiel incrémental par rapport aux champs de déplacement et d'endommagement ainsi qu'aux coordonnées nodales du maillage. Il s'agit donc d'une étude variationnelle du maillage. Nous observons que, lorsque l'endommagement atteint sa valeur maximale, l'optimisation réduit à zéro la taille de l'élément le plus endommagé. Cet élément particulier permet une représentation précise du saut de déplacement lorsque la barre se rompt. La solution globale s'avère également beaucoup plus précise qu'avec un maillage fixe. Ce travail s'inscrit dans une exploration des capacités des maillages extrêmes en mécanique numérique (X-Mesh).

**Keywords.** Phase-field, lip-field, sharp crack, X-Mesh.

**Mots-clés.** *Phase-field*, *lip-field*, fissure franche, *X-Mesh*.

**Funding.** This project has received funding from the European Research Council (ERC) under the European Union's Horizon research and innovation programme (grant agreement No. 101 071 255).

Manuscript received 18 June 2025, revised 10 July 2025, accepted 11 July 2025.

<footnote><sup>a</sup> Institute of Mechanics, Materials and Civil Engineering (iMMC), avenue Georges Lemaître 4, 1348 Louvain-la-Neuve, Belgium</footnote>

<footnote><sup>b</sup> École Centrale de Nantes, GeM Institute, UMR CNRS 6183, 1 rue de la Noë, 44321 Nantes, France</footnote>

<footnote><sup>c</sup> Nantes Université, GeM Institute, 2 rue de la Houssinière, 44322 Nantes Cedex 3, France</footnote>

<footnote>*E-mails*: nicolas.moes.remacle@uclouvain.be, benoit.le@ec-nantes.fr, nicolas.chevaugeon@univ-nantes.fr, jean-francois.remacle@uclouvain.be</footnote>

<footnote>*Corresponding author</footnote>
# 1. Introduction

The phase-field model as outlined by [1–4] has become a popular tool for studying material failure in recent years. It is based on minimising an incremental potential that depends on damage and displacement fields. One advantage of this model is that the crack path does not need to be known in advance. Cracks appear automatically during simulation and are located along the path where the damage $d$ has reached its maximum value (defined as $d = 1$ in this paper). Another similar fracture model is the lip-field model [5–7]. This model differs from the phase-field model in two ways. Firstly, the objective incremental potential does not involve the damage gradient, making it equivalent to a non-regularised potential. Secondly, a constraint is added that requires the damage field to be Lipschitz. The use of a Lipschitz constraint on the damage field can be traced back to [8]. Finally, note that the lip-field model is similar to the graded damage model introduced by [9].

The above regularised models often overlook the fact that the displacement field is discontinuous across the crack path. This can be clearly observed, for instance, in analytical one-dimensional solutions of the phase-field model. In this paper, we plan to investigate, among other things, whether this property is retained when a discrete model is considered, either with a fixed or an optimised mesh. By optimised meshes, we mean those based on the variational r-adaptation introduced in [10–21]. Variational r-adaptation is closely linked to configurational forces, as discussed in [22–26]. In the context of mesh optimisation, these forces vanish at the optimal node positions. These concepts have been applied to discontinuous fracture models, such as “Griffith-like” models [27–29] or cohesive models [30,31]. Note that mesh optimisation based strategies have also been used to improve the accuracy of phase-field results [32,33], but not in a variational framework (h-adaptativity).

Many models have attempted to incorporate discontinuities within continuous models using phase-field methods [34–37], the Thick Level Set approach [38–40] or other non-local damage models [41–45]. However, most of these models rely on specific discontinuous discretisation methods, such as X-FEM [46], and require the geometry of the macro-cracks to be identified from the continuous damage field. In this paper, we will demonstrate that displacement discontinuities can be naturally obtained by optimising the mesh using classical finite elements.

This work builds on that presented in [47–49], which used extremely deformed meshes (X-Mesh) to track moving fronts. X-Mesh is a tracking approach with a twist. The twist is that it is not always the same nodes that carry the front. Some nodes join the front, while others leave it as it evolves. The mesh can become locally highly distorted, hence the name “X-Mesh”. In the aforementioned works, X-Mesh was used for problems where the derivative of the main quantity of interest was discontinuous; for example, the temperature field in the Stefan phase-change problem [47]. Since the classical finite element method only provides approximations that are continuous at the element boundaries, discontinuities can be captured by meshing the interface.

The paper is organised as follows. Section 2 details the one-dimensional phase-field and lip-field models chosen for the study. The analytical solutions of both models are also presented. Section 3 presents the discrete models, and the finite element solutions obtained on fixed meshes are described. Section 4 describes how the mesh is optimised, highlighting the improvement in the quality of the discrete solution. This is illustrated by some numerical results in Section 5. In particular, we analyse how and when the most damaged element becomes zero-sized. The paper concludes with a discussion and some thoughts for future work in Section 6.

# 2. Phase-field and lip-field models for the fracture of a bar

We consider the deformation of a one-dimensional bar, as depicted in Figure 1 which occupies a one-dimensional domain $\Omega = [-\frac{L}{2}, +\frac{L}{2}]$. A non-decreasing bar elongation $U(t)$, rising from zero,
<figure><img src="image_3.png" /><figcaption>**Figure 1.** A bar under tension with a symmetrical damage profile.</figcaption></figure>

is imposed on the bar in a symmetric fashion:

$$u\left(\pm\frac{L}{2}, t\right) = \pm\frac{U(t)}{2}. \qquad (1)$$

The deformation of the bar at any point $x \in \Omega$ and any time $t$ is denoted by $u(x, t)$ and the corresponding strain by $\varepsilon(x, t)$. Under the small strain assumption, the strain is the derivative of the deformation: $\varepsilon = u_{,x}$. What follows considers a quasi-static evolution of $U(t)$, allowing the dependency on the time variable $t$ to be dropped and inertial effects to be neglected. Material degradation is represented by a damage variable, $d \in [0, 1]$; $d = 0$ corresponds to the undamaged material, $0 < d < 1$ to the damaged material and $d = 1$ to the fully damaged material. As material degradation is an irreversible process, $d$ cannot decrease over time. The displacement field $u$ and the damage field $d$ are solutions to a global minimisation problem

$$(u, d) = \arg \min_{u' \in \mathcal{U}, d' \in \mathcal{D}} F(u', d') \qquad (2)$$

where $F$ is the incremental potential. The specific expression of $F$ depends on whether a phase-field or lip-field model is considered. We will use the subscripts “P” and “L” to distinguish between the two. If a subscript is not present, it means that the property is valid for both. It gives the following expressions of $F_P$ and $F_L$

$$F_P(u, d) = \int_{\Omega} \frac{1}{2} \omega_P(d) E \varepsilon^2(u) \, dx + \frac{G_c}{c_P \ell_c} \int_{\Omega} \alpha_P(d) + \ell_c^2 |\nabla d|^2 \, dx, \qquad (3)$$

$$F_L(u, d) = \int_{\Omega} \frac{1}{2} \omega_L(d) E \varepsilon^2(u) \, dx + \frac{G_c}{c_L \ell_c} \int_{\Omega} \alpha_L(d) \, dx. \qquad (4)$$

The regularisation length is denoted $\ell_c$ and the scaling parameter $c$ is such that

$$c_P = 4 \int_0^1 \sqrt{\alpha_P(\beta)} \, d\beta, \quad c_L = 2 \int_0^1 \alpha_L(\beta) \, d\beta. \qquad (5)$$

The first term in the potential energy equation represents the stored elastic energy, whereas the second term represents the dissipated energy. In generic form, we will write:

$$F(u, d) = \int_{\Omega} \frac{1}{2} \omega(d) E \varepsilon^2(u) \, dx + \frac{G_c}{c \ell_c} \int_{\Omega} \alpha(d) + r \ell_c^2 |\nabla d|^2 \, dx \qquad (6)$$

where $r_P = 1$ and $r_L = 0$. The displacement field $u$ belongs to the space of the kinematically admissible solutions

$$\mathcal{U} = \{u \in H^1(\Omega) : u\left(-\frac{L}{2}\right) = -\frac{U}{2}, u\left(\frac{L}{2}\right) = \frac{U}{2}\} \qquad (7)$$

whereas

$$\mathcal{D}_P = \{d \in L^\infty(\Omega) : \underline{d} \le d \le 1\}, \qquad (8)$$

$$\mathcal{D}_L = \{d \in L^\infty(\Omega) : \underline{d} \le d \le 1 \text{ and } |d(x) - d(y)| \le |x - y|/\ell_c, \forall x, y \in \Omega\}, \qquad (9)$$
<figure><img src="image_4.png" /><figcaption>Figure 2. Phase-field and lip-field $\alpha(d)$ and $\omega(d)$ functions, for $\gamma = 0.1$.</figcaption></figure>

where $d$ is the previously computed damage value. The spaces $\mathcal{W}$ and $\mathcal{D}$ detail the conditions on the displacement and damage fields, respectively. The degradation function $\omega(d)$ is such that

$$ \omega(0) = 1, \quad \omega(1) = 0, \quad \frac{\partial \omega}{\partial d}(d) \le 0, \quad \forall d \in [0, 1] \qquad (10) $$

whereas the dissipation function $\alpha(d)$ satisfies

$$ \alpha(0) = 0, \quad \alpha(1) = 1, \quad \frac{\partial \alpha}{\partial d}(d) \ge 0, \quad \forall d \in [0; 1]. \qquad (11) $$

Finally $E$ is the Young's modulus and $G_c$ the material toughness. For the above models, the damage profile depends only on the maximum elongation reached so far. By imposing symmetry, the damage profile peaks at $x = 0$ with the peak value denoted $d_0$. Since we consider an increasing elongation starting from zero, we expect the response to first be elastic (no damage), then to show damage growth and finally to cease growing when $d_0$ reaches 1. Regarding the choice of the $\alpha$ and $\omega$ functions, we consider for $\alpha_P$ and $\omega_P$ the model detailed in [50] for the phase-field and a new model for lip-field illustrated in Figure 2:

$$ \omega_P(d) = \frac{(1-d)^2}{(1-d)^2 + 2\alpha_P(d)/(\pi\gamma)}, \quad \alpha_P(d) = 2d - d^2, \quad c_P = \pi, \qquad (12) $$

$$ \omega_L(d) = \frac{(1-d^2)^2}{(1-d^2)^2 + 2\alpha_L(d)/\gamma}, \quad \alpha_L(d) = d, \quad c_L = 1, \qquad (13) $$

where we have used

$$ \gamma = \frac{\ell_c}{\ell_{ch}}, \quad \ell_{ch} = \frac{EG_c}{\sigma_c^2}. \qquad (14) $$

A value of $\gamma$ below $\frac{8}{3\pi}$ ensures the convexity of $\omega_P$ with respect to $d$, while the maximum value of $\gamma$ for $\omega_L$ is $\frac{1}{2}$. Both are equivalent to a linear cohesive zone model (see Figure 3)

$$ \sigma = \sigma_c \left(1 - \frac{w}{w_c}\right), \quad w_c = \frac{2G_c}{\sigma_c}. \qquad (15) $$

By "equivalent", we mean that the analytical solution to minimisation problem (2) computed using the above functions or a linear cohesive zone model will produce the same elastic and dissipated energies. Usually, $\sigma_c$ and $G_c$ are material parameters that can be obtained through experiments. Then, the cohesive equivalence is valid whatever the value of $\ell_c$, provided that the value of $\gamma$ in (14) is small enough to ensure the convexity of $\omega$, which defines an upper bound on $\ell_c$.

Figure 4 shows the analytical displacement and damage fields, as well as the force versus imposed displacement curves. The detailed expressions can be found in Appendix A and B. The
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Linear cohesive law. The area under the triangle corresponds to the toughness $G_c\sigma_c w_c/2$. As the opening $w$ reaches the critical value $w_c$, the stress reaches zero.</figcaption></figure>

qualitative differences in the damage profiles between phase-field and lip-field are interesting to note (Figure 4(c)). With phase-field the width of the damaged zone is fixed, regardless of the imposed displacement $U$. In contrast, with lip-field the width of the damaged zone increases with $U$. Also, with phase-field, the damage profile is smooth at $x = 0$, except at bar breakage where it has a kink (with a slope of $\pm 1/\ell_c$). The damage $d_0$ at the center of the bar begins at a critical stress denoted by $\sigma_c$, which corresponds to a bar elongation of $U_c$. Then, damage gradually increases with elongation. Ultimately, damage reaches 1 at the critical elongation $w_c$:

$$d_0 = \begin{cases} 0 & \text{if } U \le U_c, \\ \frac{U-U_c}{w_c-U_c} & \text{if } U_c \le U \le w_c, \\ 1 & \text{if } U \ge w_c \end{cases} \qquad U_c = L \frac{\sigma_c}{E}, \qquad (16)$$

while for lip-field,

$$d_0 = \begin{cases} 0 & \text{if } U \le U_c, \\ \sqrt{\frac{U-U_c}{w_c-U_c}} & \text{if } U_c \le U \le w_c, \\ 1 & \text{if } U \ge w_c, \end{cases} \qquad U_c = L \frac{\sigma_c}{E}. \qquad (17)$$

For the models presented above to be exactly equivalent to a linear cohesive zone model, the entire damaged zone must lie strictly within the bar. At the same time, we will focus on cases with no snap-back, meaning damage must increase with an imposed increase in displacement. Combining these two conditions ((16) and (65)) gives the following bounds for the bar length for phase-field:

$$\pi \ell_c \le L \le E w_c / \sigma_c, \qquad (18)$$

which can be reformulated as

$$1 \le \frac{L}{\pi \ell_c} \le \frac{2}{\pi \gamma}, \qquad (19)$$

while for lip-field,

$$1 \le \frac{L}{\ell_c} \le \frac{2}{\gamma}. \qquad (20)$$

For longer bars, a snap-back phenomenon occurs, whereby the damage evolution is discontinuous with respect to the elongation. For shorter bars, the expression is more complex because damage growth is active over the entire length of the bar. Finally, in both models, the stress depends solely on the maximum damage value $d_0$:

$$\sigma_P = \sigma_c(1 - d_0), \qquad (21)$$

while for lip-field,

$$\sigma_L = \sigma_c(1 - d_0^2). \qquad (22)$$
<figure><img src="image_6.png" /><figcaption>(a) Non-dimensional stress vs non-dimensioned imposed displacement.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>(b) Non-dimensional displacement.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(c) Damage.</figcaption></figure>

**Figure 4.** Analytical solution. The phase-field and lip-field models produce the same global response (a) but differ terms of displacement (b) and damage profiles (c) (the curves in Figures (b) and (c) are obtained for the same imposed displacement *U*).
<figure><img src="image_9.png" /><figcaption>**Figure 5.** Problem discretisation.</figcaption></figure>

**3. The discrete model**

Figure 5 illustrates the discretisation of the problem. We consider an odd number of elements $N_e = 2n+1$. The central element has the ID 0. The nodes are located at $x_i$, $i = \pm 1, \dots, \pm(n+1)$. The set of nodal displacement and damage degrees of freedom are denoted by $u$ and $d$, respectively. The set of element sizes is denoted by $h$.

$$u = \{u_i, i = \pm 1, \dots, \pm(n+1)\}, \quad d = \{d_i, i = \pm 1, \dots, \pm(n+1)\}, \quad h = \{h_i, i = 0, \pm 1, \dots, \pm n\}. \quad (23)$$

As stated in Section 2, we will consider a symmetric evolution of damage. Therefore, we can treat the above quantities as unknown for $i = 0, \dots, n$, and then find their values for $i = -1, \dots, -n$ by symmetry. Using a one-point integration scheme, the discrete potential is given by

$$F^h(u, d) = \frac{1}{2} E \omega(d_0) \frac{(u_1 - u_{-1})^2}{h_0} + \frac{G_c}{c \ell_c} h_0 \alpha(d_0) + 2 \sum_{i=1}^{n} \left( \frac{1}{2} E \omega(\bar{d}_i) \frac{(u_{i+1} - u_i)^2}{h_i} + \frac{G_c}{c \ell_c} \left( h_i \alpha(\bar{d}_i) + r \ell_c^2 \frac{(d_{i+1} - d_i)^2}{h_i} \right) \right) \quad (24)$$

where

$$\bar{d}_i = \begin{cases} \frac{d_{-1} + d_1}{2} = d_0, & i = 0, \\ \frac{d_i + d_{i+1}}{2}, & i = 1, \dots, n. \end{cases} \quad (25)$$

The discrete potential must be minimised within the following constraints:

$$u_{n+1} = \frac{U}{2}, \quad 0 \le d_i \le 1, \quad h_i \ge 0, \quad h_0 + 2 \sum_{i=1}^{n} h_i = L, \quad i = 0, \dots, n. \quad (26)$$

By optimising with respect to the displacement field, the stress is found to be uniform with a value of

$$\sigma = E K^h(d) \frac{U}{L}, \quad K^h(d) = L \left( h_0(\omega(d_0))^{-1} + 2 \sum_{i=1}^{n} h_i (\omega(\bar{d}_i))^{-1} \right)^{-1}. \quad (27)$$

We can then eliminate the displacement field in order to obtain

$$F^h(d, U) = F_e(U) K^h(d) + G_c W^h(d), \quad (28)$$

where

$$F_e(U) = \frac{1}{2} \frac{EU^2}{L} \quad (29)$$

and

$$W^h(d) = \frac{h_0 \alpha(d_0)}{c \ell_c} + \frac{2}{c \ell_c} \sum_{i=1}^{n} \left( h_i \alpha(\bar{d}_i) + r \ell_c^2 \frac{(d_{i+1} - d_i)^2}{h_i} \right). \quad (30)$$
Regarding the constraints that $d$ and $h$ must satisfy, we define the following spaces:

$$ \mathcal{S}^h = \{(d, h) : h_0 + 2 \sum_{i=1}^{n} h_i = L\}, \qquad (31) $$

$$ \mathcal{J}^h = \{(d, h) : \underline{d}(x_i) \le d_i \le 1 \text{ and } \underline{d}_i \le d(\underline{x}_i) \le 1, \ i = 0, \dots, n\}, \qquad (32) $$

$$ \mathcal{L}^h = \{(d, h) : |d_{i+1} - d_i| \le \ell_c h_i, \ i = 1, \dots, n\}. \qquad (33) $$

The first constraint (31) imposes that the sum of the element lengths must be equal to the bar length. The second constraint (32) imposes irreversibility of the damage evolution. This amounts to requiring that the current nodal damage values $d_i$ are greater than the previous damage values interpolated at the nodes of the current mesh $\underline{d}(x_i)$, and that the current damage values interpolated at the nodes of the previous mesh $d(\underline{x}_i)$ are greater than the previous nodal damage values $\underline{d}_i$. The third constraint (33) is the Lipschitz constraint of the lip-field model. The discrete versions of the spaces $\mathcal{D}_P$ and $\mathcal{D}_L$ introduced in (8) and (9) become

$$ \mathcal{D}_P^h = \mathcal{J}^h, \qquad (34) $$

$$ \mathcal{D}_L^h = \mathcal{J}^h \cap \mathcal{L}^h. \qquad (35) $$

Finally, the discrete problem is, for a given imposed displacement $U$,

$$ (d, h) = \arg \min_{(d', h') \in \mathcal{D}^h \cap \mathcal{J}^h} F^{h'}(d', U). \qquad (36) $$

The problem is solved for $\{i = 0, 1, \dots, n\}$ using the `scipy.minimize` Python routine [51]<fnref n="1" />. Once $d$ and $h$ have been computed, $u$ can be obtained by solving an elasticity problem for which the damage is known.

## 4. X-Mesh analysis

The discrete solution of the one-dimensional fracture problem is obtained by minimising (28) with respect to the damage $d$. In this section, we will analyse what happens when minimising with respect not only to $d$, but also with respect to the mesh element size $h$. Mesh optimisation consists in imposing

$$ \frac{\partial}{\partial h_i} \left( F^h + \lambda \left( h_0 + 2 \sum_{i=1}^{n} h_i - L \right) \right) = 0, \quad i = 0, \dots, n, \qquad (37) $$

where $\lambda$ is a Lagrange multiplier that enforces the condition $h_0 + 2 \sum_{i=1}^{n} h_i = L$. This leads to

$$ -\frac{1}{2} E \left(\frac{U}{L}\right)^2 (K^h)^2 (\omega(\bar{d}_i))^{-1} + \frac{G_c}{c \ell_c} \left(\alpha(\bar{d}_i) - r \ell_c^2 \left(\frac{d_{i+1} - d_i}{h_i}\right)^2\right) + \lambda = 0. \qquad (38) $$

The Lagrange multiplier remains the same on each element. Provided at least one element remains undamaged, the value of $\lambda$ is

$$ \lambda = -(1/2) E (K^h)^2 (U/L)^2 \qquad (39) $$

In practice the length of the bar is such that the cohesive equivalence is valid, and therefore it is always possible to have an undamaged element (see end of Section 2). Furthermore, we force damage to localise on the central element using as a first damage guess in the optimisation process of the first loading step the analytical solution. This leads into

$$ \frac{1}{2} E \left(\frac{U}{L}\right)^2 (K^h)^2 \left(1 - \left(\omega(\bar{d}_i)\right)^{-1}\right) + \frac{G_c}{c \ell_c} \left(\alpha(\bar{d}_i) - r \ell_c^2 \left(\frac{d_{i+1} - d_i}{h_i}\right)^2\right) = 0. \qquad (40) $$

<footnote n="1">¹See also: https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimze.minimize.html.</footnote>
Applying the above to the central element, we have

$$ \frac{1}{2} E \left(\frac{U}{L}\right)^2 (K^h)^2 \left(1 - (\omega(d_0))^{-1}\right) + \frac{G_c}{c \ell_c} \alpha(d_0) = 0, \quad (41) $$

or, using the stress expression,

$$ \frac{1}{2} E^{-1} \sigma^2 \left(1 - (\omega(d_0))^{-1}\right) + \frac{G_c}{c \ell_c} \alpha(d_0) = 0, \quad (42) $$

giving for the phase-field model

$$ \frac{\sigma}{\sigma_c} = (1 - d_0), \quad (43) $$

while the same computation gives for the lip-field model

$$ \frac{\sigma}{\sigma_c} = (1 - d_0^2). \quad (44) $$

The X-Mesh enables the exact stress-maximum damage relationship to be recovered for both models. This is an important result. Through mesh optimisation, we can see that the stress-damage relationship corresponds to the reference relationship, regardless of the value of $d_0$ or the number of elements (provided at least one element remains undamaged). Note that if the equivalence with a linear cohesive zone model is not satisfied, for instance if we take a bar that is shorter than $2\ell_c$, we cannot use (39). Therefore all the equations up to (43) and (44) are invalid. This makes sense, since these two equations stem from the equivalence with the linear cohesive zone model.

Applying (40) to any element other than the central one gives us the following for the phase field:

$$ \frac{\ell_c^2}{h_i^2} (d_{i+1} - d_i)^2 = \alpha(\bar{d}_i) \left(1 - \left(\frac{1-d_0}{1-\bar{d}_i}\right)^2\right) = (H(\bar{d}_i, d_0))^2, \quad i=1, \dots, n-1. \quad (45) $$

The damage gradient over each element is thus related to the average damage value of the element, as defined by the reference relation (see (62) in Appendix A).

Regarding lip-field, we will now prove that for X-Mesh, the damage gradient may only be 0 or $\pm 1/\ell_c$, and not any intermediate values as it is the case for a fixed mesh. We will prove this result for any element $i > 0$, since the $i = 0$ case, which corresponds to the central element, is trivial ($d_1 = d_{-1}$ by symmetry so the gradient is necessarily zero). Let us introduce the Lagrange multiplier $\mu_i$ which imposes the Lipschitz condition

$$ \mu_i \ge 0, \quad \frac{|d_i - d_{i+1}|}{h_i} - \frac{1}{\ell_c} \le 0, \quad \mu_i \left( \frac{|d_i - d_{i+1}|}{h_i} - \frac{1}{\ell_c} \right) = 0, \quad i = 1, \dots, n, \quad (46) $$

so mesh optimisation gives

$$ \frac{\partial}{\partial h_i} \left( F^h + \lambda \left( h_0 + 2 \sum_{i=1}^n h_i - L \right) + \sum_{i=1}^n \mu_i \left( \frac{|d_i - d_{i+1}|}{h_i} - \frac{1}{\ell_c} \right) \right) = 0, \quad i = 1, \dots, n, \quad (47) $$

$$ \frac{1}{2} E \left(\frac{U}{L}\right)^2 (K^h)^2 \left(1 - (\omega(\bar{d}_i))^{-1}\right) + \frac{G_c}{\ell_c} \alpha(\bar{d}_i) - \frac{\mu_i |d_i - d_{i+1}|}{h_i^2} = 0, \quad i = 1, \dots, n, \quad (48) $$

which can be simplified using (27) and (44) to

$$ \frac{G_c}{\ell_c} \alpha(\bar{d}_i) \left(1 - \frac{(1 - d_0^2)^2}{(1 - \bar{d}_i^2)^2}\right) - \frac{\mu_i |d_i - d_{i+1}|}{h_i^2} = 0, \quad i = 1, \dots, n. \quad (49) $$

In what follows, we assume that $h_i > 0, \forall i > 0$, and demonstrate that the above equation implies that the damage gradient is either 0 or $\pm 1/\ell_c$. We will distinguish between four cases.

* $d_i = d_{i+1} = 0$: then we have that $\bar{d}_i = 0$, so $\alpha(\bar{d}_i) = 0$. As the two terms in (49) are null, this case is possible, corresponding to a fully elastic element where the gradient of $d$ is zero.
<table><caption>**Table 1.** Dimensions and material properties for numerical example.</caption><thead><tr><th>Properties</th><th>Units</th><th>Symbol</th><th>Value</th></tr></thead><tbody><tr><td>Bar length</td><td>m</td><td>L</td><td>0.2</td></tr><tr><td>Regularization length</td><td>m</td><td>l<sub>c</sub></td><td>L/5 = 0.04</td></tr><tr><td>Young's modulus</td><td>Pa</td><td>E</td><td>3 · 10<sup>10</sup></td></tr><tr><td>Fracture toughness</td><td>N/m</td><td>G<sub>c</sub></td><td>120</td></tr><tr><td>Critical tensile stress</td><td>Pa</td><td>σ<sub>c</sub></td><td>3 · 10<sup>6</sup></td></tr></tbody></table>

* $d_i = d_{i+1} = d_0$: this gives $\bar{d}_i = d_0$, so

$$1 - \frac{(1 - d_0^2)^2}{(1 - \bar{d}_i^2)^2} = 0.$$

The two terms in (49) are null which makes this case possible. This could correspond to a homogeneous evolution of damage on the bar, for instance.

* $d_i = d_{i+1} \neq 0$ and $d_i = d_{i+1} \neq d_0$: since $d_i = d_{i+1} \neq 0$, we have $\alpha(\bar{d}_i) > 0$, and $d_i = d_{i+1} \neq d_0$ gives

$$1 - \frac{(1 - d_0^2)^2}{(1 - \bar{d}_i^2)^2} \neq 0.$$

The first term of (49) is not zero, but the second term in (49) is zero, which is not possible.

* $d_i \neq d_{i+1}$, $d_i > 0$ or $d_{i+1} > 0$, $d_i < d_0$ or $d_{i+1} < d_0$. Since $d_0 > \bar{d}_i$, $\forall i > 0$, then

$$1 - \frac{(1 - d_0^2)^2}{(1 - \bar{d}_i^2)^2} > 0$$

so the first term of the above equation is positive. Since $d_i \neq d_{i+1}$, $\mu_i$ must be strictly positive for the above equation to be zero. Combining this last condition with (46) imposes that the gradient of the damage field must be $\pm 1/\ell_c$ on element $i$.

# 5. Numerical experiments

The values of the different parameters for the numerical experiments are given in Table 1. Unless otherwise specified, the number of elements per width of the fully damaged zone (which is equal to $\pi\ell_c/2$ for phase-field and $\ell_c$ for lip-field) is $n_c = 5$. Analytical curves will be plotted in black, plain lines, while numerical results will be plotted in coloured, dotted lines.

## 5.1. Phase-field analysis

Some first results obtained using phase-field are presented in Figure 6. The stress $\sigma$, computed using (27), is plotted as a function of the imposed displacement $U$ in Figure 6(a) (note that each dot corresponds to an increment of $U$). It shows that, with a fixed mesh, the bar never breaks; that is to say, the stress never reaches zero even at high values of $U$. Analysis of the discrete scheme shows that infinite elongation is required to break the bar. Conversely, with X-Mesh the bar breaks, albeit for an imposed displacement slightly smaller than the analytical value of $U = w_c$. The displacement fields are plotted in Figure 6(b). With a fixed mesh, the displacement jump when $U = w_c$ cannot be captured. However, for X-mesh, the size of the central element decreases to zero, enabling a displacement jump to be obtained. It can be seen that not only do the central nodes move, but the others do too, tending to move towards areas of higher displacement gradient. Zooming in on Figure 7 emphasises this phenomenon. Figure 6(c) shows the damage distribution on the bar. With X-Mesh, the central element has a damage value of 1,
which is not the case with a fixed mesh. This explains why the bar breaks with X-Mesh but not with a fixed mesh.

Figure 6(d) shows the $L^2$ error defined as

$$ \text{err}_2 = \sqrt{\frac{\int_{\Omega} |u - u_{\text{ex}}|^2 d\Omega}{\int_{\Omega} |u_{\text{ex}}|^2 d\Omega}} \qquad (50) $$

where $u$ is the numerical solution of the problem and $u_{\text{ex}}$ the analytical one (given by (66) or (69)). Using a middle point integration, it is computed as

$$ \text{err}_2^h = \sqrt{\frac{\sum_{i=-n}^{n} h_i |\bar{u}_i - \bar{u}_{\text{ex},i}|^2}{\sum_{i=-n}^{n} h_i |\bar{u}_{\text{ex},i}|^2}} \qquad (51) $$

where $\bar{u}_i$ and $\bar{u}_{\text{ex},i}$ are the numerical and exact displacement fields evaluated at the middle of element $i$. As can be seen, for $U > w_c$, the error obtained with X-Mesh is exactly zero, confirming the results shown in Figure 6(b). However, the error levels obtained with X-Mesh are not only lower when the bar is broken, but also for $U > w_c$. There is an error peak just before $U$ reaches $w_c$ which corresponds to the sudden stress drop shown in Figure 6(a), which we will attempt to explain later.

Figure 6(e) shows the dissipated energy as

$$ W_d(t) = \int_0^t \sigma \dot{U} dt - \frac{1}{2} \int_{\Omega} E\omega(d)\varepsilon^2 d\Omega. \qquad (52) $$

The first term in the above equation represents the energy supplied to the system, while the second term represents the stored elastic energy. The discrete counterpart is computed using a midpoint rule for spatial integration and a trapezoidal rule for temporal integration

$$ W_d^N = \sum_{k=1}^{N} \frac{\sigma^k + \sigma^{k-1}}{2} (U^k + U^{k-1}) - F_e(U^N) K^h(d^N). \qquad (53) $$

As an exception, in the above formula we noted quantities at the current time step with a $N$ exponent, and quantities at the previous time steps with a "$k = 0, \dots, N-1$" exponent.

Up to a certain imposed displacement the dissipated energy for the fixed mesh and for X-Mesh are the same. Then as $U$ gets close to $w_c$, the dissipated energy for the fixed mesh continues to increase and becomes higher than $G_c$. This is in agreement with the numerical damage profile which is wider than the analytical one shown in Figure 6(c). It also seems to keep increasing even for values of $U$ greater than $w_c$. Conversely, with X-Mesh, the dissipated energy stabilises to a value which is much closer to $G_c$.

To check the convergence of the results, the stress versus imposed displacement curves are plotted for different values of $n_c$ in Figure 8. Even with a very fine fixed mesh, the bar never breaks, whereas the results with X-Mesh are almost independent of the mesh element size. We also made a "zoom" around $U = w_c$ by adding more increments of $U$ in Figure 9. This provides a clearer view of the stress drop phenomenon observed in Figure 6(a). These figures show that as $U$ increases but remains far from $w_c$, $\sigma$ gradually decreases, following the analytical solution curve. However, at a certain value of $U$, $\sigma$ suddenly drops to zero (corresponding to $d_0$ suddenly "jumping" from a value strictly less than 1 to 1). Using smaller increments of $U$ confirms that this drop occurs abruptly and discontinuously. The evolution of $(h, u, d)$ is not continuous with respect to $U$ for X-Mesh. For some values of $U^*$, a jump occurs. As the mesh is refined, $U^*$ is closer and closer to $w_c$ and the jump gets closer and closer to zero. A possible explanation for this phenomenon will be provided in a later section.

Finally, Figure 10 shows the stress $\sigma/\sigma_c$, computed using (27) as a function of $1-d_0$. As can be seen, for X-Mesh this curve is a linear curve, confirming a property of the numerical solution that
<figure><img src="image_10.png" /><figcaption>(a) Stress versus imposed displacement (phase-field).</figcaption></figure>

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /><figcaption>(b) Displacement fields (phase-field).</figcaption></figure>

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /><figcaption>(c) Damage fields (phase-field).</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>(d) $L_2$ error on displacement (phase-field).</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>(e) Dissipated energy (phase-field).</figcaption></figure>

**Figure 6.** Phase-field: results obtained with a fixed mesh and X-Mesh, for five elements per $\frac{\pi}{2} \ell_c$.
<figure><img src="image_17.png" /><figcaption>**Figure 7.** Phase-field: nodes displacement during X-Mesh computation.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Figure 8.** Phase-field: stress versus imposed displacement curves obtained with a fixed mesh (left) and for X-Mesh (right) for different values of $n_c$.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Figure 9.** Phase-field: zoom around $U = w_c$ during X-Mesh simulations for different values of $\hat{I}$. The stress discontinuously drops to zero for $U^* < U$.</figcaption></figure>

was established in Section 4, namely that the stress satisfies the analytical relation (21) when the mesh is optimised, which is not the case with a fixed mesh.
<figure><img src="image_20.png" /><figcaption>Figure 10. Phase-field: stress $\sigma$ as a function of $d_0$.</figcaption></figure>

## 5.2. Lip-field analysis

This section presents the results obtained using the lip-field approach. The stress versus imposed displacement curves and the displacement and damage fields are plotted in Figure 11. Figure 12 shows the nodes trajectories during the X-Mesh simulation. These results are similar to those obtained using phase-field: the bar breaks with X-Mesh, but not with a fixed mesh. The size of the central elements goes to zero as $d_0$ goes to 1. The X-Mesh results in Figure 11(a) show that as with phase-field, the evolution of $(h, u, d)$ is also discontinuous. Some differences can also be observed compared to phase-field. For the fixed mesh results, the curve in Figure 11(a) shows some kind of oscillations. In fact, closer inspection of these oscillations reveals that the ascending phases correspond to elastic reloadings where damage does not evolve. This is confirmed by Figure 11(e) which shows that dissipation does not evolve during certain time steps. We will explain this phenomenon in what follows. Note that the X-mesh curve does not exhibit this issue. Also, unlike phase-field where the damage field has a zero slope at $x = 0$, except when $d_0 = 1$ when it has a kink point, the damage field with lip-field has a kink throughout the entire breaking process. (The X-Mesh recovers the features of the analytical solution mentioned in Section 2). Therefore, the mesh nodes tend to create a central element of very small size from the beginning of the simulation, reaching a zero size when $d_0 = 1$ (Figures 11(b) and (c)). With phase-field, $h_0$ decreases to zero more progressively.

Finally, with a fixed mesh, the damage slope is $\pm 1/\ell_c$ on all elements of the damaged zone (elements for which $\bar{d} > 0$), except on the outermost elements of the damage zone. We will refer to these elements as $\hat{i}$ in what follows. They are the only elements such that $\bar{d}_i > 0$, but $d$ is zero on one of their two nodes. Conversely, with X-Mesh, one node of $\hat{i}$ moves such that the slope of $d$ is exactly $\pm 1/\ell_c$, which confirms the analysis made in Section 4.

To check convergence, the stress versus imposed displacement curves are plotted for different values of $n_c$ in Figure 13. As with phase-field, the results obtained with lip-field converge to the analytical curve as the mesh is refined. In particular, with a fixed mesh, the amplitude of the elastic reloading decreases as the size of the mesh elements decreases. Also, Figure 14 shows that with X-Mesh, the stress satisfies the analytical equation (22).

To understand the phenomenon of elastic reloading observed in Figure 11(a), we focus on the behaviour of the solution around one of these elastic reloading in Figure 15. Each coloured dot in Figure 15(a) corresponds to an increment of $U$, which has a damage field counterpart of the same colour on 15(b) (which is a zoom around $\hat{i}$). It seems that the reloading (green dots in Figure 15(a)) occurs when damage becomes non-zero on a new finite element. As
<figure><img src="image_21.png" /><figcaption>(a) Stress versus imposed displacement (lip-field).</figcaption></figure>

<figure><img src="image_22.png" /></figure>

<figure><img src="image_23.png" /><figcaption>(b) Displacement fields (lip-field).</figcaption></figure>

<figure><img src="image_24.png" /></figure>

<figure><img src="image_25.png" /><figcaption>(c) Damage fields (lip-field).</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>(d) $L_2$ error on displacement (lip-field).</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>(e) Dissipated energy (lip-field).</figcaption></figure>

**Figure 11.** Lip-field: results obtained with a fixed mesh and X-Mesh, for five elements per $l_c$.
<figure><img src="image_28.png" /><figcaption>**Figure 12.** Lip-field: nodes displacement during X-Mesh computation.</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>**Figure 13.** Lip-field: stress vs imposed displacement curves obtained with a fixed mesh (left) and for X-Mesh (right) for different values of $n_c$.</figcaption></figure>

<figure><img src="image_30.png" /><figcaption>**Figure 14.** Lip-field: stress $\sigma$ as a function of $d_0$.</figcaption></figure>

expected, this corresponds to an increment where damage does not evolve (all the green curves are superimposed in Figure 15(b)). Also, during the elastic reloading, it seems that the damage slope on $\hat{i}$ is exactly $\pm\ell_c$. To gain a deeper understanding of this phenomenon, we can observe
<figure><img src="image_31.png" /><figcaption>**Figure 15.** Lip-field: elastic reloading phenomenon observed with lip-field on fixed meshes. (a) Zoom on stress vs imposed displacement curve around elastic reloading. (b) Damage field around the last element where $d > 0$. (c) Plots of $F(d_0)$. (d) Plots of $F(d_0)$, zoomed in around the discontinuity.</figcaption></figure>

that all the damage fields plotted in Figure 11(b) can be expressed as

$$d_i = \left\langle d_0 - \frac{x_i}{\ell_c} \right\rangle_+ \qquad (54)$$

where $\langle \cdot \rangle_+$ is the positive part of “■”. By substituting (54) into (28), we can express the objective function $F$ as a function of one scalar variable $d_0$ and plot it for a given value of $U$, as shown in Figures 15(c) and (d), which is a zoom around the particular value of $d_0$ that corresponds to the elastic reloading in this particular example. The function $F(d_0)$ is continuous with respect to $d_0$, but its derivative is not. It can be seen that $F(d_0)$ has the shape of a “basin”, but with some corners. As the value of $U$ increases, the “basin” rotates. For certain values of $U$ (the orange and red dots in Figures 15(c) and (d)), the minimum is on a smooth part of the basin, so its moves continuously inside its bottom. However, for certain values of $U$, the minimum becomes “stuck” in a corner of the basin, causing $d_0$ to stagnate despite $U$ increasing.

### 5.3. Lip-field five elements analysis

To understand the “jump” observed in the previous section when $d_0$ tends to 1, we will consider a simple lip-field X-Mesh example with only five elements, as illustrated in Figure 16. We assume that damage is non-zero only at the nodes of the central element. According to the analysis
<figure><img src="image_32.png" /><figcaption>**Figure 16.** Lip-field five elements case.</figcaption></figure>

made in Section 4, the damage slope must be exactly $\pm 1/\ell_c$ when the mesh is optimised, which establishes a relationship between $d_0$ and $h_1$. Additionally, the sum of the lengths of the mesh elements must equal the length of the bar, yielding the following relations:

$$h_1 = d_0 \ell_c, \quad h_2 = L/2 - (h_0/2 + d_0 \ell_c). \qquad (55)$$

The following expression for the function $F$ can be obtained for a five-element mesh, denoted by $F_5$, which depends only on $d_0$ and $h_0$:

$$F_5(d_0, h_0) = F_e(U) K_5(d_0, h_0) + \frac{G_c}{\ell_c} (h_0 \alpha_L(d_0) + 2d_0 \ell_c \alpha_L(d_0/2)), \qquad (56)$$

$$K_5(d_0, h_0) = L(h_0 \omega_L^{-1}(d_0) + 2d_0 \ell_c \omega_L^{-1}(d_0/2) + L - (h_0 + 2d_0 \ell_c))^{-1}. \qquad (57)$$

For the purposes of what follows, we will assume that $\gamma = \frac{1}{2}$. We can then express $F_5$ as a function of $d_0$ only. To achieve this, we first use the expression given in (22) and the definition of $\sigma$:

$$EU = \sigma_c (1 - d_0^2) L K_5^{-1} \qquad (58)$$

which gives (after some calculations) the expression of $h_0 d_0$:

$$h_0 d_0 = \frac{EU(1 - d_0^2)}{4\sigma_c} - \frac{L}{4}(1 - d_0^2)^2 - \frac{\ell_c d_0^2 (1 - d_0^2)^2}{(1 - (d_0/2)^2)^2} \qquad (59)$$

while the expression of $F$ simplifies to

$$F_5(d_0, h_0 d_0) = \frac{1}{2} U \sigma_c (1 - d_0^2) + 2 \frac{\sigma_c^2}{E} (h_0 d_0 + d_0^2 \ell_c). \qquad (60)$$

Combining the two equations above gives us the result that $F_5$ is a function of $d_0$ only. Assuming a linear evolution of $d$, we can determine the expression of $F_\infty$, which corresponds to the expression of $F$ that would be obtained by discretising the bar with an infinite number of elements. This is obtained by analytically integrating the different terms of $F$ in (4):

$$F_\infty(d_0) = F_e(U) \frac{L\omega(d_0)}{h_0 + 2d_0 \ell_c \frac{\omega(d_0)}{\omega(d_0/2)} + \omega(d_0)(L - (h_0 + 2d_0 \ell_c))} + G_c d_0^2. \qquad (61)$$

We plot $F_5(d_0, h_0)$, $F_5(d_0)$ and $F_\infty(d_0)$ for the values of the numerical parameters given in Table 2 in Figure 17, at five different stages as the imposed displacement $U$ increases.

(a) The first stage $U < U_c$ corresponds to no damage ($d_0 = 0$). Interestingly, Figure 17(a) shows that the corresponding $h_0$ is 0, which is sufficient to represent the behaviour of an undamaged, elastic bar.

(b) A second stage where $F_5$ has only one global minimum $d_{0,\min} < 1$. As $U$ increases, the value of $d_{0,\min}$ increases gradually and continuously.
<table><caption>**Table 2.** Dimensions and material properties for lip-field five elements example.</caption><thead><tr><td>Properties</td><td>Units</td><td>Symbol</td><td>Value</td></tr></thead><tbody><tr><td>Bar length</td><td>m</td><td>L</td><td>0.22</td></tr><tr><td>Regularization length</td><td>m</td><td>l<sub>c</sub></td><td>0.2</td></tr><tr><td>Young's modulus</td><td>Pa</td><td>E</td><td>3 · 10<sup>10</sup></td></tr><tr><td>Fracture toughness</td><td>N/m</td><td>G<sub>c</sub></td><td>120</td></tr><tr><td>Critical tensile stress</td><td>Pa</td><td>σ<sub>c</sub></td><td>3 · 10<sup>6</sup></td></tr></tbody></table>

(c) A third one where $F_5$ has two local minima: one for $d_{0,\min} < 1$ and one for $d_0 = 1$. At this stage $d_{0,\min}$ remains the global minimum.

(d) A fourth stage where $F_5$ still has two local minima $d_{0,\min} < 1$ and $1$, but now the global minimum is now at $d_0 = 1$.

(e) A fifth stage where $F_5$ only has one minimum for $d_0 = 1$.

The resolution during stage (b) will not cause any particular problems, since $F_5$ only has one global minimum. For stage (c), it should be possible to achieve a continuous increase in $d_{0,\min}$, if the solver used always finds the global minimum. However, with stage (d), $d_0 = 1$ suddenly becomes the global minimum of $F_5$. This problem can be avoided by searching for a minimum in the neighbourhood of the previous solution, thus ensuring continuity in the evolution of $d_{0,\min}$ between stages (c) and (d). The problematic stage is stage (e): indeed, the local minimum $d_{0,\min}$ “disappears”. The solution then jumps from $d_{0,\min} < 1$ to $d_0 = 1$. Note that this occurs for a value of $U$ which is still far from $w_c$. Not only does $d_0$ suddenly go to 1, but also $h_0$ suddenly goes to 0 as can be seen in Figures 17(d) and (e), left. This explains why the stress suddenly drops to zero in the results of the previous sections, assuming that this reasoning is also valid for phase-field, and for any number of finite elements in the mesh. Note that $F_\infty$ do not have this issue: whatever the value of $U$, there is only one global minimum, which increases continuously as $U$ increases.

## 6. Discussion and conclusions

In this paper, we focus on a one-dimensional fracture problem modelled using phase-field and lip-field approaches. Both approaches consist of finding the minimum of an incremental potential function with respect to displacement and damage fields. The difference lies in how the regularisation length $l_c$, which is necessary to avoid spurious mesh dependency, is introduced: through a gradient term in the expression of the incremental functional for phase-field and by a Lipschitz restriction on the damage variable for lip-field. This work examined using the X-Mesh approach to enhance the outcomes achieved with phase-field and lip-field. This approach involves moving the nodes of the finite element mesh; the movement of the nodes is determined by considering the nodal coordinates as an unknown in the optimisation problem. From this one-dimensional study, we were able to draw the following conclusions.

* Using a fixed mesh, the bar does not break for both phase-field and lip-field (i.e. the damage on the central element $d_0$ never goes to 1). Conversely, with X-Mesh, $d_0$ goes to one, while the size of the central element $h_0$ goes to zero.

* It is not only the nodes of the central element that move, but also all the elements of the mesh. This improves the quality of the solution globally, not only when $d_0 \to 1$. In particular, with the lip field, it is possible to achieve a slope in the damage field of exactly $\pm 1/l_c$, eliminating the elastic reloadings that can be observed with a fixed mesh.

* A relatively low number of elements achieves very good accuracy. In particular, the stress versus imposed displacement relationship is almost independent of the mesh element
<figure><img src="image_33.png" /><figcaption>**Figure 17.** Results for the lip-field five elements example. Left: $F_5$ as a function of $d_0$ and $h_0$ (56). Right: $F_5$ as a function of $d_0$ only (obtained using (59) and (60) and $F_\infty$. The blue square corresponds to the global minimum of $F_5$, and the black square to the global minimum of $F_\infty$.</figcaption></figure>
<figure><img src="image_34.png" /><figcaption>**Figure 17 (cont.).** Results for the lip-field five elements example. Left: $F_5$ as a function of $d_0$ and $h_0$ (56). Right: $F_5$ as a function of $d_0$ only (obtained using (59) and (60) and $F_\infty$. The blue square corresponds to the global minimum of $F_5$, and the black square to the global minimum of $F_\infty$.</figcaption></figure>

size. Using a moving mesh also allows the exact stress versus $d_0$ relation to be recovered for both phase-field and lip-field.

* However, with X-Mesh, the evolution of $d_0$ with respect to the imposed displacement $U$ is not continuous, dropping suddenly to one for a value of $U$ smaller than the critical value. This drop is due to the non-convexity of the incremental potential to optimise, which leads to the appearance/disappearance of local minima. It should be noted that this drop occurs at a value of $U$ closer to the critical value as the mesh is refined, indicating some degree of convergence towards the continuous model, for which the evolution of $d_0$ is continuous.

Finally, the main contribution of X-Mesh to fracture modelling is to adjust the finite element mesh nodes so that the material can be fully broken and displacement jumps can be modelled, which is not possible with a fixed mesh. This involves elements of size zero, a phenomenon that arises naturally from optimising the incremental potential with respect to mesh element sizes. Future work will focus on extending to 2D/3D; the main challenge here is to avoid negative Jacobian elements, which is trivial in 1D but more complex in 2D/3D. Although X-Mesh allows
for better simulation with a given number of nodes than a fixed mesh, this comes at the cost of a much slower solver. Optimising node locations is rather slow (X-Mesh computations require between two and ten times as many iterations as fixed mesh computations), and we would like to draw attention to a similar observation by Maugin in their state-of-the-art paper [25]: “One cannot avoid mentioning here problems posed by the rather slow convergence of the implied iterative procedure”. However, once a suitable solver has been found, we expect the nodes to move naturally to the location of the crack, forming layers of zero-size elements. Additionally, the stress drop phenomenon observed in 1D must be carefully studied.

# Appendix

In this section we present the analytical solutions of the phase-field and lip-field models described in Section 2.

## Appendix A. Phase-field

Regarding the damage profile over the bar, it is given implicitly by [50]

$$ \frac{dd}{dx} = -\frac{1}{\ell_c} H(d, d_0), \quad H(d, d_0) = \sqrt{(2d - d^2)\left(1 - \left(\frac{1-d_0}{1-d}\right)^2\right)}. \qquad (62) $$

The profile $d_P(x)$ is obtained by inverting the above equation. An explicit solution exists for $d_0 = 1$:

$$ d_P(x) = 1 - \sin\left(\frac{x}{\ell_c}\right). \qquad (63) $$

Also, it turns out that the width of the damaged zone obtained as

$$ \ell = \ell_c \int_{0}^{d_0} H^{-1}(\tilde{d}, d_0) d\tilde{d} \qquad (64) $$

does not depend on $d_0$:

$$ \ell = \frac{\pi}{2} \ell_c. \qquad (65) $$

This means that, as soon as the damage starts, it will occupy its full extent, which will not vary. We note that the gradient of $d$ at $x = 0$ is always 0 except when $d_0 = 1$ for which the slope is $-1/\ell_c$. Using (21) and (62), we obtain the displacement field for $U_c \le U \le w_c$ and $d_0 < 1$:

$$ u_P(x) = \begin{cases} -\frac{\sigma_c(1-d_0)}{E} \left( \ell_c \int_0^{d_0} \frac{1}{H(d, d_0)\omega(d)} dd - x - \ell \right) & \text{if } -L/2 < x < -\ell, \\ -\frac{\sigma_c \ell_c(1-d_0)}{E} \int_{d_P(x)}^{d_0} \frac{1}{H(\tilde{d}, d_0)\omega(\tilde{d})} d\tilde{d} & \text{if } -\ell < x < 0, \\ \frac{\sigma_c \ell_c(1-d_0)}{E} \int_{d_P(x)}^{d_0} \frac{1}{H(\tilde{d}, d_0)\omega(\tilde{d})} d\tilde{d} & \text{if } 0 < x < \ell, \\ \frac{\sigma_c(1-d_0)}{E} \left( \ell_c \int_0^{d_0} \frac{1}{H(d, d_0)\omega(d)} dd + x - \ell \right) & \text{if } \ell < x < L/2, \end{cases} \qquad (66) $$

while for $d_0 = 1$:

$$ u_P(x) = \begin{cases} -\frac{w_c}{2} & \text{if } x < 0, \\ \frac{w_c}{2} & \text{if } x > 0. \end{cases} \qquad (67) $$
## Appendix B. Lip-field

The damage profile is

$$d_L(x) = \begin{cases} d_0 - \frac{|x|}{\ell_c} & \text{if } |x| \le d_0 \ell_c, \\ 0 & \text{if } |x| > d_0 \ell_c, \end{cases} \qquad (68)$$

and the displacement field is for $U_c \le U \le w_c$ and $d_0 < 1$:

$$u_L(x) = \begin{cases} -\frac{\sigma_c(1-d_0^2)}{E} \left[ \ell_c \left( d_0 - d_L(x) + \frac{1}{\gamma} \left( \frac{1}{1-d_0^2} - \frac{1}{1-d_L^2(x)} \right) \right) - x - d_0 \ell_c \right] & \text{if } -L/2 < x < -d_0 \ell_c, \\ -\frac{\sigma_c(1-d_0^2)}{E} \ell_c \left[ d_0 - d_L(x) + \frac{1}{\gamma} \left( \frac{1}{1-d_0^2} - \frac{1}{1-d_L^2(x)} \right) \right] & \text{if } -\ell < x < 0, \\ \frac{\sigma_c(1-d_0^2)}{E} \ell_c \left[ d_0 - d_L(x) + \frac{1}{\gamma} \left( \frac{1}{1-d_0^2} - \frac{1}{1-d_L^2(x)} \right) \right] & \text{if } 0 < x < d_0 \ell_c, \\ \frac{\sigma_c(1-d_0^2)}{E} \left[ \ell_c \left( d_0 - d_L(x) + \frac{1}{\gamma} \left( \frac{1}{1-d_0^2} - \frac{1}{1-d_L^2(x)} \right) \right) + x - d_0 \ell_c \right] & \text{if } d_0 \ell_c < x < L/2, \end{cases} \qquad (69)$$

while for $d_0 = 1$:

$$u_L(x) = \begin{cases} -\frac{w_c}{2} & \text{if } x < 0, \\ \frac{w_c}{2} & \text{if } x > 0. \end{cases} \qquad (70)$$

## Declaration of interests

The authors do not work for, advise, own shares in, or receive funds from any organization that could benefit from this article, and have declared no affiliations other than their research organizations.

## References

[1] B. Bourdin, G. A. Francfort and J.-J. Marigo, "Numerical experiments in revisited brittle fracture", *J. Mech. Phys. Solids* **48** (2000), no. 4, pp. 797-826.

[2] A. Karma, D. Kessler and H. Levine, "Phase-field model of mode III dynamic fracture", *Phys. Rev. Lett.* **87** (2001), no. 4, article no. 045501 (4 pages).

[3] B. Bourdin, G. A. Francfort and J.-J. Marigo, "The variational approach to fracture", *J. Elasticity* **91** (2008), pp. 5-148.

[4] C. Miehe, F. Welschinger and M. Hofacker, "Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field FE implementations", *Int. J. Numer. Methods Eng.* **83** (2010), no. 10, pp. 1273-1311.

[5] N. Moës and N. Chevaugeon, "Lipschitz regularization for softening material models: the Lip-field approach", *C. R. Méc.* **349** (2021), no. 2, pp. 415-434.

[6] N. Chevaugeon and N. Moës, "Lipschitz regularization for fracture: The Lip-field approach", *Comput. Methods Appl. Mech. Eng.* **402** (2022), article no. 115644 (20 pages).

[7] N. Moës, B. Lé and A. Stershic, "Fragmentation analysis of a bar with the Lip-field approach", *Mech. Mater.* **172** (2022), article no. 104365 (20 pages).

[8] G. Dal Maso and F. Iurlano, "Fracture models as $\Gamma$-limits of damage models", *Commun. Pure Appl. Anal.* **12** (2013), no. 4, pp. 1657-1686.

[9] N. Valoroso and C. Stolz, "Graded damage in quasi-brittle solids", *Int. J. Numer. Methods Eng.* **123** (2022), no. 11, pp. 2467-2498.

[10] G. M. McNeice and P. V. Marcal, "Optimization of finite element grids based on minimum potential energy", *J. Eng. Ind.* **95** (1973), no. 1, pp. 186-190.

[11] C. A. Felippa, "Optimization of finite element grids by direct energy search", *Appl. Math. Modelling* **1** (1976), no. 2, pp. 93-96.

[12] C. A. Felippa, "Numerical experiments in finite element grid optimization by direct energy search", *Appl. Math. Modelling* **1** (1977), no. 5, pp. 239-244.
[13] E. Kuhl, H. Askes and P. Steinmann, "An ALE formulation based on spatial and material settings of continuum mechanics. Part 1: Generic hyperelastic formulation", *Comput. Methods Appl. Mech. Eng.* **193** (2004), no. 39-41, pp. 4207-4222.

[14] H. Askes, E. Kuhl and P. Steinmann, "An ALE formulation based on spatial and material settings of continuum mechanics. Part 2: Classification and applications", *Comput. Methods Appl. Mech. Eng.* **193** (2004), no. 39-41, pp. 4223-4245.

[15] J. Mosler and M. Ortiz, "On the numerical implementation of variational arbitrary Lagrangian-Eulerian (VALE) formulations", *Int. J. Numer. Methods Eng.* **67** (2006), no. 9, pp. 1272-1289.

[16] M. G. Zielonka, M. Ortiz and J. E. Marsden, "Variational r-adaption in elastodynamics", *Int. J. Numer. Methods Eng.* **74** (2008), no. 7, pp. 1162-1197.

[17] M. Scherer, R. Denzer and P. Steinmann, "On a constraint-based regularization technique for configurational r-adaptivity and 3D shape optimization", in *IUTAM Symposium on Progress in the Theory and Numerics of Configurational Mechanics* (P. Steinmann, ed.), IUTAM Bookseries, vol. 17, Springer, 2008, pp. 11-25.

[18] S. K. Lahiri, J. Bonet and J. Peraire, "A variationally consistent mesh adaptation method for triangular elements in explicit Lagrangian dynamics", *Int. J. Numer. Methods Eng.* **82** (2010), no. 9, pp. 1073-1113.

[19] J. J. Muñoz, J. P. Hambleton and S. W. Sloan, "R-adaptivity in limit analysis", in *Advances in direct methods for materials and structures* (O. Barrera, A. Cocks and A. Ponter, eds.), Springer, 2017, pp. 73-84.

[20] T. M. Tyranowski and M. Desbrun, "R-adaptive multisymplectic and variational integrators", *Mathematics 7* (2019), no. 7, article no. 642 (52 pages).

[21] Z. Shi and J. P. Hambleton, "An r-h adaptive kinematic approach for 3D limit analysis", *Comput. Geotech.* **124** (2020), article no. 103531.

[22] M. Braun, "Configurational forces induced by finite-element discretization", *Proc. Est. Acad. Sci., Phys. Math.* **46** (1997), no. 1/2, pp. 24-31.

[23] M. Braun, "Configurational forces in discrete elastic systems", *Arch. Appl. Mech.* **77** (2007), pp. 85-93.

[24] P. Steinmann, M. Scherer and R. Denzer, "Secret and joy of configurational mechanics: From foundations in continuum mechanics to applications in computational mechanics", *Z. Angew. Math. Mech.* **89** (2009), no. 8, pp. 614-630.

[25] G. A. Maugin, "Sixty years of configurational mechanics (1950-2010)", *Mech. Res. Commun.* **50** (2013), pp. 39-49.

[26] K. Schmitz and A. Ricoeur, "Theoretical and computational aspects of configurational forces in three-dimensional crack problems", *Int. J. Solids Struct.* **282** (2023), article no. 112456.

[27] C. Miehe and E. Gürses, "A robust algorithm for configurational-force-driven brittle crack propagation with R-adaptive mesh alignment", *Int. J. Numer. Methods Eng.* **72** (2007), no. 2, pp. 127-155.

[28] M. Scherer, R. Denzer and P. Steinmann, "Energy-based r-adaptivity: A solution strategy and applications to fracture mechanics", *Int. J. Fract.* **147** (2007), pp. 117-132.

[29] A. Qinami, E. C. Bryant, W. C. Sun and M. Kaliske, "Circumventing mesh bias by r- and h-adaptive techniques for variational eigenfracture", *Int. J. Fract.* **220** (2019), pp. 129-142.

[30] G. Geißler, C. Netzker and M. Kaliske, "Discrete crack path prediction by an adaptive cohesive crack model", *Eng. Fract. Mech.* **77** (2010), no. 18, pp. 3541-3557.

[31] U. De Maio, F. Greco, P. Lonetti and A. Pranno, "A combined ALE-cohesive fracture approach for the arbitrary crack growth analysis", *Eng. Fract. Mech.* **301** (2024), article no. 109996 (17 pages).

[32] Y. Li, W. Lai and Y. Shen, "Variational h-adaption method for the phase field approach to fracture", *Int. J. Fract.* **217** (2019), pp. 83-103.

[33] F. Freddi and L. Mingazzi, "Adaptive mesh refinement for the phase field method: A FEniCS implementation", *Appl. Eng. Sci.* **14** (2023), article no. 100127 (12 pages).

[34] B. Giovanardi, A. Scotti and L. Formaggia, "A hybrid XFEM-Phase field (Xfield) method for crack propagation in brittle elastic materials", *Comput. Methods Appl. Mech. Eng.* **320** (2017), pp. 396-420.

[35] R. J. M. Geelen, Y. Liu, J. E. Dolbow and A. Rodríguez-Ferran, "An optimization-based phase-field method for continuous-discontinuous crack propagation", *Int. J. Numer. Methods Eng.* **116** (2018), no. 1, pp. 1-20.

[36] A. Muixí, O. Marco, A. Rodríguez-Ferran and S. Fernández-Méndez, "A combined XFEM phase-field computational model for crack growth without remeshing", *Comput. Mech.* **67** (2021), pp. 231-249.

[37] J. Zhang, H. Yu, W. Xu, C. Lv, M. Micheal, F. Shi and H. Wu, "A hybrid numerical approach for hydraulic fracturing in a naturally fractured formation combining the XFEM and phase-field model", *Eng. Fract. Mech.* **271** (2022), article no. 108621.

[38] N. Moës, C. Stolz, P.-É. Bernard and N. Chevaugeon, "A level set based model for damage growth: The thick level set approach", *Int. J. Numer. Methods Eng.* **86** (2011), no. 3, pp. 358-380.

[39] P.-É. Bernard, N. Moës and N. Chevaugeon, "Damage growth modeling using the Thick Level Set (TLS) approach: Efficient discretization for quasi-static loadings", *Comput. Methods Appl. Mech. Eng.* **233-236** (2012), pp. 11-27.

[40] L. A. T. Mororó, A. Poot and F. P. van der Meer, "Skeleton curve and phantom node method for the Thick Level Set approach to fracture", *Eng. Fract. Mech.* **268** (2022), article no. 108443 (24 pages).
[41] M. R. R. Seabra, P. Šuštarič, J. M. A. Cesar de Sa and T. Rodič, "Damage driven crack initiation and propagation in ductile metals using XFEM", *Comput. Mech.* **52** (2013), pp. 161–179.

[42] E. Tamayo-Mas and A. Rodríguez-Ferran, "A new continuous-discontinuous damage model: cohesive cracks via an accurate energy-transfer process", *Theor. Appl. Fract. Mech.* **69** (2014), pp. 90–101.

[43] Y. Wang and H. Waisman, "From diffuse damage to sharp cohesive cracks: A coupled XFEM framework for failure analysis of quasi-brittle materials", *Comput. Methods Appl. Mech. Eng.* **299** (2016), pp. 57–89.

[44] S. Sarkar, I. V. Singh and B. K. Mishra, "A simplified continuous-discontinuous approach to fracture based on decoupled localizing gradient damage method", *Comput. Methods Appl. Mech. Eng.* **383** (2021), article no. 113893.

[45] A. Negi and S. Kumar, "A continuous-discontinuous localizing gradient damage framework for failure analysis of quasi-brittle materials", *Comput. Methods Appl. Mech. Eng.* **390** (2022), article no. 114434.

[46] N. Moës, J. E. Dolbow and T. Belytschko, "A finite element method for crack growth without remeshing", *Int. J. Numer. Methods Eng.* **46** (1999), no. 1, pp. 131–150.

[47] N. Moës, J.-F. Remacle, J. Lambrechts, B. Lé and N. Chevaugeon, "The eXtreme Mesh deformation approach (X-MESH) for the Stefan phase change model", *J. Comput. Phys.* **477** (2023), article no. 111878 (30 pages).

[48] A. Quiriny, J. Lambrechts, N. Moës and J.-F. Remacle, "X-Mesh: A new approach for the simulation of two-phase flow with sharp interface", *J. Comput. Phys.* **501** (2024), article no. 112775.

[49] A. Chemin, J. Lambrechts, N. Moës and J.-F. Remacle, "Solving the porous medium equation with the eXtreme Mesh deformation approach (X-Mesh)", 2025. Online at https://arxiv.org/abs/2501.03083.

[50] J.-Y. Wu, "A unified phase-field theory for the mechanics of damage and quasi-brittle failure", *J. Mech. Phys. Solids* **103** (2017), pp. 72–99.

[51] P. Virtanen et al., "SciPy 1.0: fundamental algorithms for scientific computing in Python", *Nat. Methods* **17** (2020), no. 3, pp. 261–272.