See discussions, stats, and author profiles for this publication at: https://www.researchgate.

net/publication/369898266

Reconstruct lower-dimensional crack paths from phase-field point cloud

Article  in  International Journal for Numerical Methods in Engineering · April 2023

DOI: 10.1002/NME.7249

CITATION READS

1 285

3 authors:

Yue Xu Tao You

Hohai University Helmholtz-Zentrum für Umweltforschung
1 PUBLICATION   1 CITATION    12 PUBLICATIONS   96 CITATIONS   

SEE PROFILE SEE PROFILE

Qizhi Zhu
Hohai University
135 PUBLICATIONS   2,972 CITATIONS   

SEE PROFILE

Some of the authors of this publication are also working on these related projects:

Contitutive Modelling of Geomaterials View project

Micromechanics-enhance phase-field model for quasi-brittle materials View project

All content following this page was uploaded by Tao You on 08 April 2023.

The user has requested enhancement of the downloaded file.

Received <day> <Month>, <year>; Revised <day> <Month>, <year>; Accepted <day> <Month>, <year>

DOI: xxx/xxxx

ARTICLE TYPE

Reconstruct lower-dimensional crack paths from phase-field point

cloud

Yue Xu1,3 | Tao You*1,2,3 | Qizhi Zhu1,3

1 Key Laboratory of Geomechanics and

Embankment Engineering of Ministry of Abstract
Education, Hohai University, Nanjing
210098, China
Because of the smeared representation of phase-field fracture, reconstructing lower-
2 Department of Environmental Informatics, dimensional crack paths is challenging, which may impede the further applications
Helmholtz Centre for Environmental of the phase-field method in the modeling of fault friction and fluid flow in frac-
Research – UFZ, Leipzig 04318, Germany
3 College of Civil and Transportation ture. In the present work, we propose to capture the crack curves or surfaces from
Engineering, Hohai University, Nanjing two-dimensional (2D) or three-dimensional (3D) phase-field point cloud using an
210098, China
optimized ridge regression algorithm, and k-nearest neighbor (kNN) and principal
Correspondence component analysis (PCA) are used to estimate the normal direction of each segment
*Corresponding author name, Corresponding
on the identified discrete crack path. The sensitivity and computational efficiency
address. Email: tao.you@ufz.de
of the proposed method are investigated thoroughly. Subsequently, this method is
extended to reconstruct the complex discrete fracture networks (DFNs). Finally, sev-
eral benchmarks for fracture, fault friction and fluid flow problems are presented to
demonstrate the strength of the proposed approach.

KEYWORDS:
Phase-field method; Crack reconstruction; Fault friction; Point cloud; Normal direction; Discrete fracture
network

1 INTRODUCTION

The precise prediction of fractures and coupling processes in fractures are of great significance for ensuring safety, produc-
tivity and durability during the exploitation of geo-energy reservoirs. Continuum damage method, as a common mathematical
approach, especially the non-local damage [35, 34], has been widely used to describe the fracture phenomena characterized
by a fracture process zone (FPZ)[41, 6]. This method produces a smeared representation of a sharp crack, which is favorable
for numerical implementation and incorporation of complex constitutive relationships. As a special gradient damage model
[14, 25], the phase-field method has become a powerful tool for modeling fracture and coupling problems with multiphysics,
e.g., thermal[29, 4], hydraulic[22], and chemical[39, 43] couplings. Nevertheless, such a smeared approximation of fracture
gives rise to the difficulty in reconstructing a lower-dimensional crack path and complex discrete fracture network (DFN), which
may hinder the further application of the phase-field model for fluid flow and fault friction.
Taking the hydraulic phase-field model as an example, the fluid flow in the fracture is distinct from that in the reservoir domain
where the former may follow lubricant law or Navier-Stokes equations, as considered in [37, 16, 11, 12]. As a result, the sharp
crack topology that is the solution domain of the fluid flow equation should be explicitly presented. Meanwhile, the normal
direction of this sharp crack is required for calculating the fracture aperture [28, 26, 45]. When considering the fault friction, on
the other hand, the sharp crack topology is also required to enable accurate calculation of the normal and tangential stresses on
the crack surface for specific friction criteria[18, 10].
2

Several hybrid phase-field models have been proposed to circumvent the above-mentioned issues in reconstructing sharp crack
topology, for example by hybridizing with extended finite element method (XFEM) [20, 47] and embedded discrete fracture
method (EDFM) [12]. These hybrid models naturally guarantee the sharp crack topology, but also inherit the disadvantage of
these two discrete models (XFEM, EDFM), i.e., it is difficult to handle complex fracture networks in 3D . On the other hand,
the level-set method [32] is an attractive approach to identifying the crack domain and estimating normal vector on the level-set
boundary[26, 45, 42]. This method can provide a set of iso-surface describing the boundary of the crack domain, but cannot
identify the accurate lower-dimensional crack path. A point searching algorithm is still needed, for example in [18] . In addition,
the normal direction estimated by the level-set method is not feasible at the crack tip because of a blurred representation by the
phase-field model, and therefore requires ad hoc numerical treatment, as illustrated by [45].
Regarding directly reconstructing the sharp crack topology in the continuum finite element simulation, Bottoni et al. [9] pro-
posed a crack-tracking method in the post-processing stage by maximizing the damage variable in a given direction at a given
distance. Tamayo-Mas and Rodríguez-Ferran [40] used a medial-axis-based model to define the crack path in the gradient-
enhanced damage model. In the phase-field model, Ziaei-Rad et al. [48] and Rudy et al. [20] presented the optimization
algorithms to capture the path within the phase-field framework in two-dimensional conditions. Fan and Choo [18] proposed a
searching algorithm to identify the discrete crack path and corresponding normal directions in order for frictional slip simula-
tion. Recently, Yang et al. [44] proposed an identifying algorithm considering the phase-field value as the indicator and cutting
physical patches to get explicit cracks. Zeng et al. [46] developed a new identifying method for multi-branched crack tips which
searches paths in local regions and updates the regions automatically.
However, these crack reconstruction methods are based on point searching algorithms which are hard to tackle the crack band
problem, i.e., the material is fully damaged in a band [5, 30, 14]. Additionally, most of the above-mentioned approaches have not
been extended to tackle three-dimensional fracture problems due to the difficulties or even limitations of the tracking algorithm
in 3D. As a result, we propose a set of numerical approaches to identify discrete crack path (i.e., curve path for 2D fracture and
surface path for 3D fracture) and the corresponding normal direction of each segment on this crack path. A thorough evaluation
of the computational efficiency, robustness and bias will be provided. The proposed method has no dimension limitation and is
well suited to reconstructing the DFNs in real time. Furthermore, the mesh-caused bias will be greatly alleviated.
The remainder of this paper is organized as follows. In Section 2, the basic formulations of the phase-field fracture model
are given first, and then the proposed methods for identification of single crack path and complex DFN follow. The estimation
approach of the normal direction is introduced. Section 3 gives several fitting results for investigating the effectiveness, sensitivity
and computational efficiency of the proposed optimization method. In Section 4, several numerical examples are presented in
which the proposed method will be used to address fracture, friction and fluid flow problems. The newly proposed method
will be compared with other existing methods in 2D fracture and friction problems (in view that 3D problems have not been
considered for other existing identification methods). The capability and advantages of the proposed method will be illustrated
with two typical benchmarks. Finally, the paper ends with a brief summary and prospect for future work.

2 METHODS

We start with the formulations of the phase-field method for rock fracture and fault friction. Following that, an optimization
method for identifying a lower-dimensional crack path from the phase-field point cloud is proposed, and the corresponding
normal direction will be calculated by principal component analysis (PCA). Finally, the optimization method will be extended
from a single discrete crack to a complex fracture network by giving a set of classifying and clustering algorithms.

2.1 Phase-field formulation for rock fracture and fault friction

Consider a cracked domain Ω ⊂ ℝ𝑛 with boundary 𝜕Ω and dimension 𝑛 = 2, 3. As shown in Fig. 1 , 𝒖𝐷 ∶ 𝜕Ω𝑢 → ℝ𝑛
and 𝒕𝑁 ∶ 𝜕Ω𝑡 → ℝ𝑛 are respectively the Dirichlet and Neumann boundary conditions which satisfy 𝜕Ω𝑢 ∪ 𝜕Ω𝑡 = 𝜕Ω and
𝜕Ω𝑢 ∩ 𝜕Ω𝑡 = ∅. The embedded crack is regularized by the well-known phase-field method which gives the crack surface density
function as [27] ( )
1 𝑑2
𝛾(𝑑, ∇𝑑) = + 𝓁|∇𝑑|2 (1)
2 𝓁
where 𝑑 is the phase-field variable ranging from 0 (intact material) to 1 (fully damaged material), and 𝓁 is the regularized length
parameter.
 3

FIGURE 1 Stress decomposition in the interface-oriented system.

2.1.1 Balance laws

The linear momentum balance law reads without considering body force

∇ ⋅ 𝝈 = 𝟎 in Ω, (2)

and the boundary conditions are as follows

{
𝒖 = 𝒖𝐷 on 𝜕Ω𝑢
(3)
𝝈 ⋅ 𝒏𝑡 = 𝒕𝑁 on 𝜕Ω𝑡
On the other hand, we employ AT2 phase-field model to capture the induced cracks, and the volumetric–deviatoric strain
decomposition [1] is used to avoid crack penetration under compression. As a result, we have the phase-field equation as follows
( )
𝑑
𝑔 ′ (𝑑) + 𝐺𝑐 − 𝓁∇ ⋅ ∇𝑑 = 0 in Ω (4)
𝓁
where 𝐺𝑐 is the critical energy release rate and  is the crack driving force which takes the form
( )
𝑘0
(𝒙, 𝑡) = max 𝜓+ [𝜺(𝒙, 𝑠)] = max ⟨tr 𝜺⟩2+ + 𝜇0 𝜺dev ∶ 𝜺dev (5)
𝑠∈[0,𝑡] 𝑠∈[0,𝑡] 2
where the bracket operator is defined as ⟨⋅⟩+ = (⋅ + | ⋅ |) ∕2 and the deviatoric strain is defined as 𝜺dev = 𝕂 ∶ 𝜺 with the
fourth order projection tensor 𝕂 (with its component defined as 𝐾𝑖𝑗𝑘𝑙 = [𝛿𝑖𝑘 𝛿𝑗𝑙 + 𝛿𝑖𝑙 𝛿𝑗𝑘 ]∕2 − 𝛿𝑖𝑗 𝛿𝑘𝑙 ∕3) for deviatoric strain and
𝜺 = 1∕2(∇𝒖 + ∇𝒖T ).

2.1.2 For fracture problem

We first focus on the fracture problem. Given the volumetric-deviatoric strain decomposition above, the degraded stress is given
as
( )
𝝈 = 𝑔(𝑑) 𝑘0 ⟨trε⟩+ 𝐈 + 2μ0 εdev + 𝑘0 ⟨trε⟩− 𝐈 (6)

where the bracket operator is defined as ⟨⋅⟩− = (⋅ − | ⋅ |) ∕2 and 𝐈 is the second-order identity tensor. 𝑘0 and 𝜇0 are the bulk
modulus and shear modulus of the material respectively. Additionally, 𝑔(𝑑) is the damage degradation function which takes the
form
𝑔(𝑑) = (1 − 𝜅)(1 − 𝑑)2 + 𝜅 (7)

in which 𝜅 is a constant and 0 < 𝜅 << 1.

Accordingly, the stress-strain tangent tensor in this case can be written as
𝜕𝝈
ℂ= = 𝐻[tr𝜺]𝑔(𝑑)𝑘0 I ⊗ I + 2𝑔(𝑑)𝜇0 𝕂 + 𝐻[−tr𝜺]𝑘0 I ⊗ I (8)
𝜕𝜺
where 𝐻(⋅) is the Heaviside function.
4

2.1.3 For friction problem

When an existing crack (fault) is closed, friction may occur and the stick-slip behavior on the crack surface should be taken into
account. To this end, we employ the method proposed by Fan and Choo [18] for frictional contact. As shown in Fig. 1 , the
Cauchy stress is decomposed into three scenarios according to the open-closed states of the fault, i.e.,
⎧ 𝑔(𝑑)𝝈 no contact
⎪ bulk
𝝈 = ⎨ 𝝈bulk stick (9)
( ( ) )
⎪ 𝝈bulk + 𝜇𝑝N, bulk sign 𝜏bulk − 𝜏bulk (𝒏 ⊗ 𝒕 + 𝒕 ⊗ 𝒏) slip
⎩
where 𝜇 denotes the friction coefficient of the interface, and the bulk stress 𝝈bulk is given as

𝝈bulk = ℂbulk ∶ 𝜺, (10)

with ℂbulk = 3𝑘0 𝕁 + 2𝜇0 𝕂 being the stiffness tensor of the material. Additionally, 𝑝N, bulk and 𝜏bulk are respectively the normal
and tangential contact stresses based on the interface-oriented coordinate system in Fig. 1 , which take the forms
{
𝑝N, bulk = −𝝈bulk ∶ (𝒏 ⊗ 𝒏)
𝜏bulk = 𝝈bulk ∶ (𝒏 ⊗ 𝒕)
Accordingly, the stress-strain tangent tensor can be derived as
⎧ 𝑔(𝑑)ℂ no contact
⎪ bulk
ℂ = ⎨ ℂbulk stick (11)
( )
⎪ ℂbulk + [1 − 𝑔(𝑑)] ℂ𝑓 − ℂ𝜏 slip
⎩
where
{ ( )[ ]
ℂ𝑓 = −𝜇 sign 𝜏bulk (𝑘0 − 23 𝜇0 )(𝒏 ⊗ 𝒕 + 𝒕 ⊗ 𝒏) ⊗ 𝐈 + 2𝜇0 (𝒏 ⊗ 𝒕 + 𝒕 ⊗ 𝒏) ⊗ (𝒏 ⊗ 𝒏)
ℂ𝜏 = 𝜇0 (𝒏 ⊗ 𝒕 + 𝒕 ⊗ 𝒏) ⊗ (𝒏 ⊗ 𝒕 + 𝒕 ⊗ 𝒏)
It turns out that the normal direction of the existing crack is crucial for tackling the friction problem. One of the common
∇𝑑
methods is using the gradient of phase-field, i.e., 𝒏 = ‖∇𝑑‖ . However, the accuracy of this approximation is insufficient due to (1)
a blunt approximation at crack tips and (2) an undesirable approximation when 𝑑 is very close to 1. In this paper, after identifying
the discrete crack path, the normal direction is readily estimated using principal component analysis and the above-mentioned
problems will be circumvented.
The numerical implementation of the above two-field (𝒖−𝑑) problem represented by Eq. (2) and (4) can be found in Appendix
A. It is noteworthy that when the volumetric-deviatoric strain decomposition is used to simulate brittle shear fracture under
compression i.e., trε < 0, the volumetric locking phenomenon will lead to instabilities when the fracture localizes, as has been
pointed out by [19]. During the numerical test of 2D single-edge notched shear fracture presented in Section 4.1, this instability
behavior will abort the computation when the fracture approaches the edge of the model. As a result, a mixed finite element
discretization is employed in this paper to satisfy the LBB condition[2, 3].

2.2 Reconstruction of single discrete crack from phase-field point cloud

In this subsection, we propose to identify the single discrete fracture from the phase-field point cloud via interpolation and
optimization techniques. After that, the normal directions for the lower-dimensional identified path will be estimated.
We start with the definition of the phase-field point cloud. As shown in Fig. 2 , the smeared crack is represented by the phase-
field variable which is defined on the mesh nodes. A set of points (i.e., point cloud) can be extracted from the mesh nodes given
a specified threshold of the phase-field variable.
Given the determined point cloud as sample data, the realistic sharp crack should be situated in the middle of these points
ideally, as shown in the rightmost of Fig. 2 . Although it is straightforward to use an interpolation technique as present in
[48, 18], a searching algorithm is needed and its application to 3D fracture is still challenging. Rather than by interpolation, we
first use an optimization strategy (inspired by RegularizeNd [31] ) to find a single lower-dimensional crack, i.e., a curve in 2D
or a surface in 3D, and then it will be employed in the reconstruction of the discrete fracture network.
 5

Volume Point cloud Crack surface

Optimization algorithm

FIGURE 2 Definition of phase-field point cloud and discrete crack path.

2.2.1 Approximation of lower-dimensional crack path

After obtaining the phase-field point cloud, we proceed to approximate the realistic sharp crack path. To explain this process,
we suppose the number of sample points (i.e., the determined point cloud) is 𝑛 and the number of output points (i.e., the nodes
on the discrete crack curve or surface) is 𝑝. In most situations the number of samples 𝑛 falls behind the number of outputs, i.e.,
𝑛 < 𝑝. The coordinates of sample data and output points are designated as (𝐗, 𝐘) and (𝜶, 𝜷) respectively. For a 2D problem, 𝐗
and 𝜶 have one column which delineates the coordinate in one of the dimensions. For the 3D problem, 𝐗 and 𝜶 will have two
columns that correspond to the coordinates in two of the dimensions.
The coordinates of output points 𝜶 are given first as a mesh grid which is spanned by 𝑝1 × 𝑝2 (= 𝑝) mesh points. As a result,
one can customize the number and density of points in the desired crack path 1 . The next step is to calculate the coordinate 𝜷,
thereby identifying the spatial position of the sharp crack. A direct approach is by interpolation, i.e.,
{
𝐗 = 𝐀𝜶
(12)
𝐘 = 𝐀𝜷
where 𝐀 is (𝑛 × 𝑝)-dimensional interpolation matrix. Taking linear interpolation in 2D as an example, the coordinate 𝐗𝑖 (𝑖 =
1, ..., 𝑛) of an arbitrary sample point can be expressed as a linear combination of its two adjacent output points 𝜶𝑗 and 𝜶𝑘
(𝑗, 𝑘 = 1, ..., 𝑝), i.e., 𝐗𝑖 = 𝑎𝜶𝑗 + 𝑏𝜶𝑘 , where the weights 𝑎 and 𝑏 satisfy 𝑎 + 𝑏 = 1. Therefore, the interpolation matrix 𝐀 can
be determined by collecting the weights for all the sample points. In 3D, a bilinear interpolation is needed and the principle is
identical to that in 2D.
Our goal here is to find the linear relation 𝐘𝑖 = 𝐀𝑖,∗ 𝜷∗ from the samples by using regression analysis. The loss function using
conventional sum-of-squares can be written as
(𝜷) = ‖𝐘 − 𝐀𝜷‖22 (13)
Let derivative of  equal to zero, one can obtain the so-called normal equation, i.e.,
𝐀T 𝐀𝜷 = 𝐀T 𝐘 (14)
The estimation of 𝜷 then reads
( )−1 T
𝜷̂ = 𝐀T 𝐀 𝐀 𝐘, (15)
( )−1
if 𝐀 𝐀 T
is well-defined.
Unfortunately, in most cases, the interpolation matrix 𝐀 is high-dimensional in view of 𝑛 < 𝑝, which makes the estimator
incalculable due to the super-collinearity [15]. A typical threat for regression analysis in this situation is overfitting as all the
covariates are modeled which brings the noise. It is also the reason why point search algorithms are needed for the other methods,
for example in [48, 18]. Alternatively, another method to circumvent super-collinearity is using a ridge regression estimator as
shown below.

1 A coarse or fine mesh grid will be a potential factor that influences the quality of the approximation of a realistic sharp crack. Therefore, we will denote ℎ
̂ min as the
minimum space size of the mesh grid and investigate its effect in Appendix B
6

2.2.2 Ridge regression estimation

The principle of ridge regression can be considered as introducing an additional penalty term in the loss function (13) of the
original linear regression problem. In this paper, we want to guarantee the smoothness of the discrete crack paths, and therefore,
define the ridge loss function as
(𝜷) = ‖𝐘 − 𝐀𝜷‖22 + 𝜆 ‖𝐃𝜷‖22 (16)
where 𝐃 is the second-derivative operator which can be readily calculated by finite difference techniques given the prescribed
mesh grid 𝜶. The operator 𝐃 is a (𝑚 × 𝑝)-dimensional matrix with 𝑚 being the number of second-derivatives for 𝑝 output points,
and we have 𝑚 = 𝑝 − 2 in 2D and 𝑚 = 𝑝1 (𝑝2 − 2) + 𝑝2 (𝑝1 − 2) in 3D. Therefore, the size of matrix 𝐃 relies on the number of
output points we predefined.
Additionally, the parameter 𝜆 in Eq. (16) is the penalty parameter which dictates the contribution of ridge penalty term to
the loss function and eventually affects the location of the discrete crack. The larger 𝜆 generates the larger contribution of the
smoothness constraint to the loss function. In particular, a large enough 𝜆 will give a line in 2D or a plane in 3D. To specify 𝜆,
we use the form similar to [31], i.e.,
𝑛
𝜆 = 𝑘𝑠
𝑚
where 𝑘𝑠 is an independent parameter to regulate the smoothness of the expected fitting path. The introduction of a scaling
term 𝑚𝑛 is to balance the residuals caused respectively by the liner regression part and penalty part (in view that the number of
residuals for the penalty part is 𝑚), thereby obtaining an objective result with varying output points.
The normal function of the regression problem in Eq. (16) takes the form
( T )
𝐀 𝐀 + 𝜆𝐃T 𝐃 𝜷 = 𝐀T 𝐘 (17)
The ridge regression estimation of 𝜷 then arrives at
( )−1 T
𝜷̂ = 𝐀T 𝐀 + 𝜆𝐃T 𝐃 𝐀 𝐘 (18)
( T )
In this estimator, a non-zero value of 𝜆 will give a positive definite matrix for 𝐀 𝐀 + 𝜆𝐃T 𝐃 [21] so as to avoid the unsolvability
due to the above-mentioned super-collinearity.
Without loss of generality, we take the one-dimensional crack path approximation as an example to reflect the advantage of
the optimized ridge regression. The comparison of identification results with and without ridge optimization is shown in Fig.
3 . The black points corresponding to the phase-field point cloud are extracted from the red damage band given the damage
threshold 𝑑c = 0.98. The green crack path using the ridge regression estimator is smoother than that using the linear regression
estimator, and it is situated in the middle of the damage band more ideally. The crack path estimated by ridge regression is more
mesh-independent, which gives a relatively objective crack geometry and location.

(a) Linear regression estimator Eq.(15) (b) Ridge regression estimator Eq.(18)

FIGURE 3 Comparison of discrete crack paths obtained from linear regression and ridge regression estimators.
 7

2.2.3 Normal direction estimation

With the discrete crack path at hand, we can estimate the local normal direction on this crack path. A normal vector 𝒏 will be
assigned to each point on the crack path by finding a hyperplane point by point. The k Nearest Neighbours (kNN)[13] is used
first to find its neighbor points and then the Principal Component Analysis (PCA)[23] will be used to determine the local plane
and the corresponding normal vector.
To concretize this process, considering a point 𝒙 on the crack path, we would like to select its k nearest neighbors in 𝑟-ball
{ }
which satisfy 𝒙𝑖 ∣ ‖ ‖
‖𝒙𝑖 − 𝒙‖ < 𝑟 , and then find a local plane Π which minimizes the sum of square distances. As such, the
problem boils down to
∑k
( )2
 ∶= min dist 𝒙𝑖 , Π (19)
𝑖=1
The distances above are real orthogonal distances between the points and the corresponding local plane.
∑k
Next, we define the centroid 𝒙𝑐 as 𝒙𝑐 = 1k 𝑖=1 𝒙𝑖 . The vector from 𝒙𝑖 to this centroid can be expressed as 𝒚𝑖 = 𝒙𝑖 − 𝒙𝑐 .
Accordingly, the original problem (19) can be recast into
∑
k
 ∶= min (𝒚𝑖𝑇 ⋅ 𝒏)2
‖𝒏‖=1
𝑖=1
( )
where 𝒏 is the unit normal vector of the plane Π. By denoting the set 𝒚1 𝒚2 ⋯ 𝒚k as 𝒀 , the problem can be further simplified
as
 ∶= min
𝑇
𝒏𝑇 (𝒀 𝒀 𝑇 )𝒏 (20)
𝒏 𝒏=1
The optimization problem is further converted into PCA of a covariance matrix regarding to the nearest neighbors for a point.
PCA gives a fitting hyperplane and a set of orthogonal bases to best describe a given data set. This algorithm is implemented
based on the eigenvalue decomposition of the covariance matrix. If we use 𝑺 to substitute 𝒀 𝒀 𝑇 , Eq. (20) turns into a constrained
minimization problem, which can be solved with Lagrange Multiplier Methods [7]
{ ( )
(𝒏, 𝜆) = 𝒏𝑇 𝑺𝒏 − 𝜆 𝒏𝑇 𝒏 − 1
∇ = 0

The equation 𝑺𝒏 = 𝜆𝒏 can be derived from the condition 𝜕

𝜕𝒏
= 0. The hyperplane normal 𝒏 is the eigenvector of 𝑺 with the
smallest eigenvalue. The matrix 𝑺 can be decomposed as
⎛𝜆1 ⎞
𝑺 = 𝑽 ⎜ ⋯ ⎟𝑽 𝑇 (21)
⎜ ⎟
⎝ 𝜆𝑑 ⎠
with 𝜆𝑖 and 𝑽 representing the eigenvalues and eigenvectors of the matrix 𝑺 respectively.
The matrix 𝑺 measures the variance of the data points along the direction 𝒗 in Fig. 4 . The big eigenvalue means the corre-
sponding eigenvector predicts the direction where the points possess strong variance. Therefore, the eigenvector related to the
smallest eigenvalue is the targeted normal direction.

Neighbor points Big eigenvalue Small eigenvalue

FIGURE 4 The correlation between data variance and eigenvalue.

8

2.3 Reconstruction of the complex discrete fracture network

In this section, we focus on reconstructing the complex DFNs. A set of algorithms is proposed to formulate the workflow in
identifying intersected crack networks and discretizing them into elements. The aforementioned algorithms in Section 2.2 will
serve in determining fracture fragments, and the new algorithm works for clustering and intersecting of fractures.

2.3.1 Classification of fracture fragments and joints

We first introduce the concept of surface variation [33] to determine whether a point in the phase-field point cloud belongs to a
fracture fragment or joint. The surface variation 𝑠𝑟 at point 𝒙 is defined as
𝜆
𝑠𝑟 (𝒙) = ∑𝑑1 (22)
𝜆
1 𝑖
where 𝜆𝑖 corresponds to the eigenvalues of the matrix 𝑺 as expressed in Eq. (21), and 𝑟 represents the range search radius for PCA
calculation. The surface variation reflects the spatial distribution of the neighboring points within the range 𝑟. The maximum
value of 𝑠𝑟 is 1/3, and the minimum value is 0.
Under the 3D condition, as shown in Fig. 5 , the maximum value of 𝑠𝑟 is 1/3 which means the neighboring points distribute
isotropically and without a preferable direction. On the contrary, the points appear to be located on a quasi-plane with a preferable
direction when the value of 𝑠𝑟 is close to 0.

FIGURE 5 Points distribution under different surface variations.

As a result, according to the value of surface variation 𝑠𝑟 , the fracture joint can be distinguished from the fracture fragments
based on the phase-field point cloud. For example, in Fig. 6 , a Y-shaped phase-field point cloud is presented and the value of 𝑠𝑟
is higher in the intersection area. By giving a threshold 𝑠𝑐𝑟 = 0.11, we can divide the phase-field point cloud into two categories,
i.e., points belonging to fracture fragments and joints.

FIGURE 6 Classifying and clustering of Y-shaped phase-field point cloud.

 9

2.3.2 Clustring for fracture fragments and joints

Since the algorithm proposed in Section 2.2 is valid only for a single discrete fracture, a clustering algorithm can then be applied
to group each fracture joint and fragment for individual treatment. In this paper, we employ the DBSCAN clustering algorithm
[17, 36, 38]. This algorithm is based on a density-based notion of clusters and designed to discover clusters of arbitrary shape.
The key idea is that for each point of a cluster, the neighborhood within a given radius has to contain at least a minimum number
of points, i.e. the density in the neighborhood has to exceed some threshold.
Given the classified phase-field point clouds for fracture fragments or fracture joints, Algorithm 1 is performed for each
unlabeled point so as to determine the core points of each cluster. Taking the Y-shaped phase-field point cloud as an example
again, in Fig. 6 , three fracture fragments can be distinguished.

Algorithm 1 for DBSCAN clustering

Input: The searching distance: 𝐸𝑝𝑠; the minimum number of neighbor points: 𝑀𝑖𝑛𝑃 𝑡𝑠.
1: Calculate the number of neighbor points of an arbitrary unlabeled point within a given distance 𝐸𝑝𝑠. If the number of the
neighbor points is less than 𝑀𝑖𝑛𝑃 𝑡𝑠, this point is treated as noise. On the contrary, this point is marked as the core point of
this cluster and assigned a new cluster label.
2: Access one of the obtained neighbor points. If it has not been assigned a cluster label, the newly created cluster label is
assigned to it.
3: Evaluate whether this neighbor point is a core point by following the method in step one. If it is the core point, step two is
repeated.
4: Repeat the second and third steps until no more core points can be found.
Output: Labeled phase-field point cloud for each fracture fragment and joint.

Subsequently, the algorithms proposed in Section 2.2 will be used to reconstruct each discrete fracture fragment. In the
intersection area, the discrete crack paths will be generated to bridge the three fracture fragments. Note that a twist 3D fracture
has not been considered in the present work, which means the crack surface is linear in its depth direction. Accordingly, the
discrete fracture paths which bridge the fracture fragments can be simplified by several intersecting planes.
Finally, we summarize the entire workflow for reconstructing complex discrete fracture networks as in Algorithm 2.

Algorithm 2 for reconstructing complex discrete fracture networks from the phase-field point cloud.
Input: Coordinates and phase-field values of finite element nodes;
1: Extract the point clouds whose phase-field value are greater than the damage threshold (e.g. 0.99) and denote them as set
𝑝𝑐 ;
2: Calculate the surface variation of each point based on Eq. (22);
3: Divide the point clouds into two categories, i.e., fracture fragments 𝑓 𝑟𝑎𝑔𝑠 and fracture joints 𝑗𝑜𝑖𝑛𝑡𝑠 , according to the value
of surface variation;
4: Perform DBSCAN clustering algorithm[17, 36] to distinguish different point cloud clusters in set 𝑓 𝑟𝑎𝑔𝑠 and set 𝑗𝑜𝑖𝑛𝑡𝑠 ;
5: Reconstruct the lower-dimensional paths for each fracture fragment and joint;
6: Estimate the normal direction of the discrete crack path using the method presented in Section 2.2.3;
Output: Discrete fracture networks and normal vectors;

3 RESULTS

In this section, we present two numerical examples to demonstrate the feasibility of the proposed workflow. The first example
corresponds to a single discrete crack in 2D and 3D respectively. In the second example, a 3D discrete fracture network is
addressed.
10

3.1 Single discrete crack in 2D and 3D

We first consider two self-designed wave-typed crack models as shown in Fig. 7 (a) and Fig. 8 (a). For the 2D case, the size
of the square specimen is 1mm × 1mm which is discretized by 30,010 plane quadrilateral elements with a minimum mesh size
ℎmin = 0.0035 mm in the region of local damage. The length scale parameter is taken as 𝓁𝑠 = 0.02 mm . For the 3D case, the size
of the specimen is 1mm × 1mm × 0.2mm which is discretized by 356,540 8-node linear hexahedron elements. The parameters
are the same as those in the 2D case.

(a) Mesh and prescribed realistic crack curve (b) Phase-field (c) Point cloud and fitting crack curve

FIGURE 7 Reconstruction of wave-typed discrete crack in 2D.

(a) Mesh and prescribed realistic crack surface (b) Phase-field (c) Point cloud and fitting crack surface

FIGURE 8 Reconstruction of wave-typed discrete crack in 3D.

First, the realistic crack path, i.e., 2D curve in Fig. 7 (a) or 3D surface in Fig. 8 (a), is prescribed in the middle of the model
by constraining 𝑑 = 1 on the points of the realistic crack path. For a linear crack, our optimization method gives an almost
identical crack path compared with the realistic one. Therefore, we increase the curvature of the crack path to demonstrate the
robustness of our optimization method. The phase-field equation is then solved to give a smeared representation of the realistic
crack, as shown in Fig. 7 (b) and Fig. 8 (b). The phase-field point cloud is shown as red scattered paths in Fig. 7 (c) and Fig.
8 (c). The lower-dimensional crack paths represented by the green curve under a prescribed grid density are approximated by
using the proposed workflow.
 11

It can be found that the discrete cracks are smooth even though the phase-field point cloud is unevenly distributed. Additionally,
the precision and efficiency of the proposed method are examined in Appendix B. The proposed workflow can reconstruct
discrete cracks in 2D and 3D conditions with satisfactory precision and few computational efforts.

3.2 Complex discrete fracture network in 3D

A complex DFN is considered in this section. The size of the specimen is 100mm × 100mm × 50mm which is discretized by
8-node hexahedron elements with a minimum mesh size ℎmin = 1 mm in the region of local damage. The length scale parameter
is taken as 𝓁𝑠 = 2 mm. The two parameters used in DBSCAN algorithm are 𝐸𝑝𝑠 = 2 mm and 𝑀𝑖𝑛𝑃 𝑡𝑠 = 8. Similar to the case
in Section 3.1, the same procedure is applied to generate phase-field crack. As shown in Fig. 9 , the deflection, intersection and
merging of fracture fragments can be correctly captured. The normal directions of the DFN are also estimated, as shown in Fig.
9 (d).

(a) Mesh and prescribed discrete fracture network (b) Phase-field value of the volume

(c) Lower-dimensional crack paths (d) Lower-dimensional crack paths and normals

FIGURE 9 Reconstruction of a complex discrete fracture network in 3D.

4 APPLICATIONS
12

4.1 Real-time reconstruction of sharp crack during fracturing

In this section, we first incorporate the above proposed workflow into the phase-field modeling of fracturing. The sharp 2D and
3D cracks will be reconstructed in real-time during crack propagation.
The single-edge notched shear test (SENS) is considered here. The geometry and boundaries of the numerical model are
presented in Fig. 10 (a). For the 3D case, the thickness of the model is taken as 0.2 mm. A horizontal displacement is applied
on the top edge with a constant loading rate of 𝑑𝑢 = 1 × 10−4 mm per loading step in the first 100 steps and 𝑑𝑢 = 5 × 10−5
mm per loading step in the next 100 steps. In the 2D case, the model is discretized by 26,059 quadratic 8-node quadrilateral
serendipity elements with a minimum size ℎmin = 0.00375 mm for local refinement. In the 3D case, the model is discretized
by 211,190 8- node hexahedron elements with the same size of refined elements. The parameters used in the simulation are as
follows: Young’s modulus 𝐸0 = 2.1 × 105 MPa, Poisson’s rate 𝑣0 = 0.3, fracture energy parameter 𝐺𝑐 = 2.7 N/mm, the length
scale parameter 𝓁𝑠 = 0.0075 mm. Note that a strict global convergence by multiple 𝑢 − 𝑑 iterations is not constrained for the
3D case and 3-5 iterations are enough to give plausible results.

(a) Geometry and boundaries (b) Step 120

(c) Step 195 (d) Step 200

FIGURE 10 The evolution of lower-dimensional crack during the fracture of 2D single-edge notched specimen. (𝑑𝑐 =
0.95, ℎ̂ min = 0.2ℎmin )

We extract the phase field at different loading steps in Fig. 10 and Fig. 11 for 2D and 3D cases respectively. The lower-
dimensional crack paths shown by the green curves in Fig. 10 are generated using the proposed reconstructing algorithm 2. The
damage threshold for extracting phase-field point cloud is taken as 𝑑𝑐 = 0.95, and the mesh grid size for sharp cracks is taken
as ℎ̂ min = 0.2ℎmin . As expected, the method proposed in this paper can effectively reconstruct the low-dimensional fracture in
real time during the fracture process.

4.2 Fault friction

In this section, we investigate the fault friction problems in which the discrete crack path and the normal direction of each point
on the crack are required for calculating the contact stresses. The first numerical case involves both fault friction and fracturing,
while the second case considers a pure fault friction problem. We adopt the same phase-field formulation for the frictional
 13

(a) Step 65 (b) Step 95 (c) Step 200

FIGURE 11 The evolution of lower-dimensional crack during the fracture of 3D single-edge notched specimen.(𝑑𝑐 =
0.95, ℎ̂ min = 0.2ℎmin )

problem as proposed by Fan and Choo [18], as expressed in Section 2.1.3. The only difference lies in how to reconstruct the
lower-dimensional crack path and estimate its normal directions.

4.2.1 Propagation of an inclined frictional crack

The geometry of the model is shown in Fig. 12 (a) where an inclined crack is embedded starting from (0.0,0.7) m to (1.3,2.0)
m. The pre-existing crack is modeled by applying a Dirichlet boundary as 𝑑 = 1. A vertical displacement loading is applied on
the top edge with a constant loading rate of 𝑑𝑢 = 2 × 10−4 m per loading step until a through-going fracture is generated. The
model is discretized by 16,100 8-node structured quadratic quadrilateral elements with a minimum size ℎmin = 0.02 mm for local
refinement along the possible crack path. The parameters used in the simulation are as follows: Young’s modulus 𝐸0 = 1 × 104
MPa, Poisson’s rate 𝑣0 = 0.3, fracture energy parameter 𝐺𝑐 = 50 N/mm, length scale parameter 𝓁𝑠 = 0.04 mm.
We extract the phase field at three different loading steps in Fig. 12 . The damage threshold 𝑑𝑐 and mesh grid size ℎ̂ min are
chosen as 0.99 and 0.2ℎmin for reconstructing the lower-dimensional crack. The discrete crack paths are depicted in Fig. 12 by
the green curves. The crack initiates from the tip of the inclined frictional crack and propagates to the right edge of the model.

(a) Geometry and boundary con- (b) Step 57 (c) Step 58 (d) Step 59
ditions

FIGURE 12 The evolution of lower-dimensional crack during loading of frictional fracture test.

Although the fracture pattern does not show much difference compared with that in [18], the resulting sharp crack from the
proposed method is smoother and the normal directions are more continuously distributed. As shown in Fig. 13 , we apply
both the piece-wise approximation [18] and our proposed method to reconstruct the discrete crack curves and the corresponding
14

normal directions under the same mesh. Clearly, the piece-wise linear approximation method is heavily influenced by the mesh,
while the crack path obtained by the proposed method in Fig. 13 (b) is almost unaffected by the mesh. This characteristic is
more significant for the pure friction problem because the prediction of the normal direction directly affects the calculation of
the contact stress on the fracture, which gives rise to different mechanical responses.

(a) The piece-wise linear approximation

(b) The proposed optimization method

FIGURE 13 Comparison of the discrete crack paths reconstructed by the piece-wise linear approximation [18] and the proposed
method. (𝑑𝑐 = 0.99, ℎ̂ min = 0.2ℎmin )

4.2.2 Friction of internal crack

We then illustrate the effect of the lower-dimensional crack paths obtained by the two different approximation methods on the
mechanical response of friction. The geometry of the model is shown in Fig. 14 where an internal crack is embedded starting
from (0.3,0.3) m to (0.7,0.7) m. The crack is not allowed to propagate in this case so that a pure friction problem is considered.
A vertical displacement loading is applied on the top edge with a constant loading rate of 𝑑𝑢 = 0.01 m per loading step.
The model is discretized by unstructured quadrilateral elements with a minimum size ℎmin = 0.004 m for local refinement. The
parameters used in the simulation are as follows: Young’s modulus 𝐸0 = 1 × 104 MPa, Poisson’s rate 𝑣0 = 0.3, fracture energy
parameter 𝐺𝑐 = 50 N/mm, length scale parameter 𝓁𝑠 = 0.016 m.
As shown in Fig. 15 , the discrete fractures are reconstructed with the piece-wise linear approximation method [18] and the
proposed method respectively. As expected, the proposed method gives a smoother crack path under this unstructured mesh.
The relative horizontal and vertical displacements (𝑢𝑥 and 𝑢𝑦 ) across the fault are presented in Fig. 16 where the values of 𝑢𝑥
and 𝑢𝑦 are normalized by their maximum values. It is clear that due to the mesh-induced roughness of the discrete crack path,
the piece-wise linear approximation causes a disturbance of the displacement field, whereas the proposed method produces a
more regular one.
 15

(a) Geometry and boundary conditions (b) Initial phase-field value

FIGURE 14 Setup of the friction problem for internal crack.

(a) The piece-wise linear approximation

(b) The proposed optimization method

FIGURE 15 Comparison of discrete fractures reconstructed by the piece-wise linear approximation method [18] and the pro-
posed method respectively. (𝑑𝑐 = 0.99, ℎ̂ min = 0.2ℎmin )

4.3 Fluid flow in discrete fracture network

In this section, we consider the fluid flow in the fractures in the phase-field model. Due to the unknown nature of discrete
fractures in the smeared damage model, most of the existing phase-field models consider an enhancement of permeability in the
fracture domain and neglect the solution of fluid flow equation on the lower-dimensional DFNs, see the review article by [22]
and the references therein. One beneficial attempt is the work by Santill et al. [37] which considers the fluid flow on a sharp
crack line. However, their approach is difficult to extend to complex DFNs.
With the proposed method above, the DFNs can be reconstructed from the phase-field model, and then the fluid flow on the
DFNs can be easily captured, for example, by using dfnWorks [24]. To do so, we consider a phase-field model as shown in Fig.
17 (a). The size of the model is 4 m × 2 m × 2 m. After reconstruction, the DFN in Fig. 17 (b) comprises a center-located
horizontal fracture with a length of 2 m and four inclined fractures dipping at an angle of 45°. The aperture of the fractures is
taken as 1 cm. The inlets and outlets are marked in Fig. 17 (c). The pressures on the inlet and outlet boundaries are set as 2 MPa
16

(a) The relative displacement in 𝑥 direction (b) The relative displacement in 𝑦 direction

FIGURE 16 The relative displacements across the fault under two different discrete fractures reconstructed by the piece-wise
linear approximation method [18] and the proposed method respectively.

and 1 MPa respectively. The flow process is considered as the single-phase saturated incompressible flow (Richards equation
[24]) with permeability and fluid viscosity being 1 × 10−12 m2 and 1 × 10−3 Pa⋅s, respectively.

(a) The phase-field model (b) Discrete fracture network

(c) Layout of inlet and outlet boundaries (d) Pressure filed

FIGURE 17 Procedures of modeling fluid flow in discrete fracture network in the phase-field model.

Fig. 17 (d) presents the resulting pressure field. The fracture flow problem in the phase-field model can be directly solved
after reconstructing the DFNs using the proposed method above. Benefiting from the accurate reconstruction of DFNs, this
method will be soon extended to tackle the hydro-mechanical coupled processes between the discrete fracture network and the
solid matrix.

5 CONCLUSION

In this work, we present a workflow for reconstructing low-dimensional fracture paths and estimating fracture normal direc-
tions based on the concept of phase-field point cloud. We provide several numerical examples in 2D and 3D to illustrate the
 17

effectiveness of the proposed method. Several advantages are highlighted, such as (1) the proposed method is feasible for 2D
and 3D complex fracture networks with arbitrary structural or nonstructural finite element meshes, (2) the normal direction of
the phase-field crack can be correctly estimated, and (3) the lower-dimensional fracture paths can be reconstructed in real time
with high computational efficiency. The proposed method can contribute to addressing the challenge of phase-field modeling
for fault friction and fluid flow in DFNs. In our future work, a hydro-mechanical coupling model that explicitly represents the
discrete fracture network will be developed to simulate the coupled processes in fractures and reservoirs, and the DFNs will be
considered to be updated in real time during the evolution of the phase field.

CREDIT AUTHOR STATEMENT

Yue Xu: Investigation, Methodology, Software, Validation, Visualization, Writing - Original Draft. Tao You: Conceptualiza-
tion, Methodology, Investigation, Software, Visualization, Writing - Original Draft, Writing - Review & Editing. Qi-zhi Zhu:
Supervision, Resources, Investigation

ACKNOWLEDGMENTS

The author Tao You would like to acknowledge the supports by the Helmholtz-OCPC Postdoc-Fellowship (No. 202) and the
National Natural Science Foundation of China (No. 12202137).

CONFLICT OF INTEREST

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to
influence that work reported in this paper.

APPENDIX

A MIXED FINITE ELEMENT DISCRETIZATION

The above 𝒖−𝑑 problem Eqs. (2) and (4) will be solved by finite element method. The weak form can be obtained by multiplying
the strong forms with wight function 𝒘𝑢 and 𝑤𝑑 , respectively, and then integrating the equations over the domain Ω, i.e.,

∇𝒘𝑛 𝝈dΩ − 𝒕𝑁 ⋅ 𝒘𝑢 dΓ = 0 (A1)

∫ ∫
Ω Γ𝑁
( )
𝑑𝑤𝑑
𝑔 ′ (𝑑)𝑤𝑑 dΩ + 𝐺𝑐 + 𝓁∇𝑑 ⋅ ∇𝑤𝑑 dΩ = 0 (A2)
∫ ∫ 𝓁
Ω Ω

with the trail solution spaces 𝑢 for displacement and 𝑑 for the phase-field defined as
{ }
𝑢 ∶= 𝒖 ∣ 𝒖 ∈ 𝐻 1 , 𝒖 = 𝒖𝐷 on 𝜕Ω𝑢 (A3)
{ }
𝑑 ∶= 𝑑 ∣ 𝑑 ∈ 𝐻 1 (A4)
where 𝐻 1 means the first order Sobolev space.
The displacement 𝒖 and phase-field 𝑑 can be approximated by Galerkin’s method as follows
𝒖 = 𝑵𝑢 𝒖,
̂ 𝒘𝑢 = 𝑵𝑢 𝒘̂ 𝑢 , ∇𝒖 = 𝑩𝑢 𝒖,
̂ ∇𝒘𝑢 = 𝑩𝑢 𝒘̂ 𝑢 (A5)
̂ 𝑤𝑑 = 𝑵𝑑 𝑤̂ 𝑑 , ∇𝑑 = 𝑩𝑑 𝑑,
𝑑 = 𝑵𝑑 𝑑, ̂ ∇𝑤𝑑 = 𝑩𝑑 𝑤̂ 𝑑 (A6)
18

Phase‐field degree of freedom
Displacement degree of freedom
Gauss quadrature point

FIGURE A1 Schematic illustration of 𝒖 − 𝑑 mixed finite element.

where the hatted variable denotes the value at nodes. 𝑵𝑢 and 𝑵𝑑 are the shape functions for displacement and phase-field sepa-
rately, and 𝑩𝑢 and 𝑩𝑑 are the corresponding derivatives of shape functions. In order to safisfy LBB conditions, the displacement
shape function is taken to correspond to the quadratic 8-node quadrilateral serendipity element (Q8), while the phase-field shape
function is taken as the bilinear 4-node quadrilateral element (C4). The schematic illustration of the mixed element formulation
is shown in Fig. A1 .
Accordingly, the displacement and the phase-field problems can be converted into the following linear system as

𝑲𝑢 𝒖̂ = 𝑭𝑢 (A7)
𝑲𝑑 𝑑̂ = 𝑭𝑑 (A8)

with

𝑲𝑢 = 𝑩𝑢T ℂ𝑩𝑢 dΩ (A9)

∫
Ω

𝑭𝑢 = 𝑵𝑢T 𝒃 dΩ + 𝑵𝑢T 𝒕𝑁 dΩ (A10)

∫ ∫
Ω 𝜕Ω𝑡

where the stress-strain tangent tensor ℂ has been defined in Eq. (8) for fracture problem and Eq .(11) for friction problem.
Additionally, for the phase-field problem, we have
{( ) }
𝐺𝑐
𝑲𝑑 = 2 + T T
𝑵𝑑 𝑵𝑑 + 𝐺𝑐 𝓁𝑩𝑑 𝑩𝑑 dΩ (A11)
∫ 𝓁
Ω

𝑭𝑑 = 2𝑵𝑑T  dΩ (A12)
∫
Ω

After that, we use the staggered algorithm to solve Eqs. A7 and A8, namely the displacement 𝒖𝑖+1 is solved first with fixed
phase-field 𝑑 𝑖 and then the phase-field 𝑑 𝑖+1 is solved with fixed displacement 𝒖𝑖+1 until the global convergence criterion is
satisfied.

B ESTIMATION OF ERRORS AND COMPUTATIONAL EFFICIENCY

In this appendix, the precision and efficiency of the proposed method for reconstructing 2D and 3D discrete cracks will be
evaluated. Regarding the examples in Section 3, two parameters, i.e., damage threshold 𝑑𝑐 and minimum element size ℎ̂ min ,
are investigated to figure out how they affect the fitting result for the realistic crack path. For the quantitative evaluation of the
 19

results, a similarity measure is proposed, i.e.,

⎧ √ ⎫
√ 𝑛 ( )2
⎪ √∑ ⎪
 ∶= ⎨min √ 𝑘 𝑘
𝑎𝑖 − 𝑏𝑗 , 𝑎𝑖 ∈ 𝐴⎬ (B13)
𝑏𝑗 ∈𝐵
⎪ 𝑘=1 ⎪
⎩ ⎭
where 𝑛 represents the dimension of the model. 𝑎𝑖 is an arbitrary point belonging to the fitting path set 𝐴 and 𝑏𝑗 is an arbitrary
point belonging to the realistic path set 𝐵. The mechanism is as follows: For each member in set 𝐵, calculate the Euclidean
distance to all the members in set 𝐴, find out the minimum of the distances which represents the closest pair, and then keep
looping until all the members in set 𝐵 have found their closest pairs. The distance for each member in set 𝐵 is defined as the
deviation distance, which is depicted in Fig. B2 .

FIGURE B2 Crack path deviation from the realistic sharp crack path with large curvature.

To describe the deviation magnitude more objectively, we define a deviation ratio as


 ∶= (B14)
ℎmin
After the aforementioned modeling and optimizing processes, we compare the realistic crack path with the fitted one by
varying the values of damage threshold 𝑑𝑐 and the minimum size of mesh grid ℎ̂ min . To clarify the deviation ratio  for different
influencing factors and keep the data as complete as possible, we employ the violin plots.
• For the 2D case:
As shown in Fig. B3 , the damage threshold ranges from 0.8 to 0.95 with an interval of 0.01. The mesh grid size ℎ̂ min varies
from 0.25ℎmin , 0.5ℎmin , 0.75ℎmin to 1ℎmin . Note that, in practice, we generally expect a fine mesh grid for the fitting crack path,
i.e., letting ℎ̂ min < ℎmin .
The deviation ratio for each closest pair between the realistic path and the fitting path is presented on the left half of the violin
column, and the right half of the violin represents the distribution of these scattered data on the left. Additionally, the white
dots represent the mean values of the deviation ratios. From the whole picture of Fig. B3 , the mean deviation distance will not
exceed half of the minimum element length. As 𝑑c increases, the deviation ratio rises first until 𝑑𝑐 = 0.85 and then falls. When
𝑑c > 0.85, as the threshold increases, the number of points in the phase-field point cloud gradually decreases and the fitting path
converges to the true path. When 𝑑c < 0.85, as the threshold decreases, the number of points on both sides of the crack becomes
more even and the fitting path tends to be the true path. Note that, a slight influence of 𝑑𝑐 is due to the smeared pattern of the
phase field, i.e., in the areas with large curvature of the crack curve, more points inside the curve are captured than the outside
given a specific value of 𝑑𝑐 . Hence the deviation ratio shows an increase in that situation. In practice, we generally expect a
relatively big value of 𝑑𝑐 which will further improve computational efficiency as verified below.
Meanwhile, with the increase of the mesh grid size, the mean deviation ratio increases gradually, indicating that the effect of
mesh grid spacing is more prominent than the damage threshold. Therefore, we need to set the value of ℎ̂ min as small as possible,
and then the computation cost should be taken into account. As shown in Fig. B4 , we present the CPU time and memory
for the crack identification processes under different values of 𝑑𝑐 and ℎ̂ min ∕ℎmin . The proposed algorithms are implemented by
20

(a) ℎ̂ min /ℎmin =0.25 (b) ℎ̂ min /ℎmin =0.5

(c) ℎ̂ min /ℎmin =0.75 (d) ℎ̂ min /ℎmin =1

FIGURE B3 Deviation ratios with different damage threshold 𝑑c and mesh grid size ℎ̂ min of the discrete crack path for the 2D
model.

Julia language [8] and BenchmarkTools.jl (https://github.com/JuliaCI/BenchmarkTools.jl) is applied to record running time and
memory.

(a) Time consumption (b) Memory occupation

FIGURE B4 Computational costs with different damage threshold 𝑑𝑐 and mesh grid size ℎ̂ min of discrete crack path in the 2D
case.
 21

In Fig. 4 (a), the time consumption rises slightly with 𝑑c decreasing from 0.98 to 0.8. Given the same damage threshold, the
computation time increases gradually as the grid spacing gets smaller. In this case, the longest time consumption for the finest
mesh grid of the fitting crack path (i.e., ℎ̂ min ∕ℎmin = 0.1) is 2 milliseconds, which is almost imperceptible comparing to the
phase-field fracture modeling. In Fig. 4 (b), the memory occupation shows similar behaviors as the time consumption when
varying 𝑑𝑐 and ℎ̂ min ∕ℎmin .

• For the 3D case:

Similar to the 2D case, we investigate the deviation ratio by varying the damage threshold and mesh grid. The range of ℎ̂ min
remains the same as that in 2D, while the damage threshold ranges from 0.7 to 0.9 with an interval of 0.02. The deviation
ratios under different damage thresholds and mesh grid sizes are presented in Fig. B5 where the deviation ratio shows a similar
variation with the 2D case and the mean value is under 50%.

(a) ℎ̂ min /ℎmin =0.25 (b) ℎ̂ min /ℎmin =0.5

(c) ℎ̂ min /ℎmin =0.75 (d) ℎ̂ min /ℎmin =1

FIGURE B5 Deviation ratios with different damage threshold 𝑑c and mesh grid size ℎ̂ min of the discrete crack path for the 3D
model.

Regarding the computational efficiency in the 3D case, as shown in Fig. 6 (a), the time slightly increases with 𝑑c decreasing
from 0.9 to 0.7, and the calculation time increases more rapidly when the mesh grid size becomes smaller. For this 3D case,
the maximum time consumption is 1.5 seconds which is also acceptable in the phase-field modeling of 3D fracture. In Figure.
6 (b), the memory occupation shows as the time consumption when varying 𝑑𝑐 and ℎ̂ min ∕ℎmin .
22

(a) Time consumption (b) Memory occupation

FIGURE B6 Computational costs with different damage threshold 𝑑𝑐 and mesh grid size ℎ̂ min of discrete crack path in the 3D
case

References

[1] Amor, H., J.-J. Marigo, and C. Maurini, 2009: Regularized formulation of the variational brittle fracture with uni-
lateral contact: Numerical experiments. Journal of the Mechanics and Physics of Solids, 57, no. 8, 1209–1229,
doi:https://doi.org/10.1016/j.jmps.2009.04.011.
URL https://www.sciencedirect.com/science/article/pii/S0022509609000659

[2] Babuška, I., 1971: Error-bounds for finite element method. Numerische Mathematik, 16, no. 4, 322–333,
doi:https://doi.org/10.1007/BF02165003.

[3] Babuška, I. and R. Narasimhan, 1997: The babuška-brezzi condition and the patch test: an example. Computer Methods
in Applied Mechanics and Engineering, 140, no. 1, 183–199, doi:https://doi.org/10.1016/S0045-7825(96)01058-4.
URL https://www.sciencedirect.com/science/article/pii/S0045782596010584

[4] Badnava, H., M. A. Msekh, E. Etemadi, and T. Rabczuk, 2018: An h-adaptive thermo-mechanical phase field model for
fracture. Finite Elements in Analysis and Design, 138, 31–47, doi:https://doi.org/10.1016/j.finel.2017.09.003.
URL https://www.sciencedirect.com/science/article/pii/S0168874X17303347

[5] Bažant, Z. P. and B. H. Oh, 1983: Crack band theory for fracture of concrete. Matériaux et construction, 16, no. 3, 155–
177, doi:https://doi.org/10.1007/BF02486267.
URL https://link.springer.com/article/10.1007/BF02486267

[6] Bazant, Z. P. and J. Planas, 2019: Fracture and size effect in concrete and other quasibrittle materials. Routledge.

[7] BERTSEKAS, D. P., 1982: Chapter 2 - the method of multipliers for equality constrained problems. Constrained Opti-
mization and Lagrange Multiplier Methods, D. P. BERTSEKAS, ed., Academic Press, 95–157.
URL https://www.sciencedirect.com/science/article/pii/B9780120934805500064

[8] Bezanson, J., A. Edelman, S. Karpinski, and V. B. Shah, 2017: Julia: A fresh approach to numerical computing. SIAM
review, 59, no. 1, 65–98, doi:https://doi.org/10.1137/141000671.

[9] Bottoni, M., F. Dufour, and C. Giry, 2015: Topological search of the crack pattern from a continuum mechanical compu-
tation. Engineering Structures, 99, 346–359, doi:https://doi.org/10.1016/j.engstruct.2015.05.005.
URL https://www.sciencedirect.com/science/article/pii/S0141029615003168

[10] Bryant, E. C. and W. Sun, 2021: Phase field modeling of frictional slip with slip weakening/strengthening under non-
isothermal conditions. Computer Methods in Applied Mechanics and Engineering, 375, 113557.
 23

[11] Chukwudozie, C., B. Bourdin, and K. Yoshioka, 2019: A variational phase-field model for hydraulic
fracturing in porous media. Computer Methods in Applied Mechanics and Engineering, 347, 957–982,
doi:https://doi.org/10.1016/j.cma.2018.12.037.
URL https://www.sciencedirect.com/science/article/pii/S0045782518306443

[12] Costa, A., M. Cusini, T. Jin, R. Settgast, and J. E. Dolbow, 2022: A multi-resolution approach to hydraulic fracture
simulation. International Journal of Fracture, 237, no. 1-2, 165–188.

[13] Cover, T. and P. Hart, 1967: Nearest neighbor pattern classification. IEEE Transactions on Information Theory, 13, no. 1,
21–27, doi:10.1109/TIT.1967.1053964.

[14] de Borst, R. and C. V. Verhoosel, 2016: Gradient damage vs phase-field approaches for fracture:
Similarities and differences. Computer Methods in Applied Mechanics and Engineering, 312, 78–94,
doi:https://doi.org/10.1016/j.cma.2016.05.015, phase Field Approaches to Fracture.
URL https://www.sciencedirect.com/science/article/pii/S0045782516303796

[15] Draper, N. R. and H. Smith, 1998: Applied regression analysis, volume 326. John Wiley & Sons.

[16] Ehlers, W. and C. Luo, 2018: A phase-field approach embedded in the theory of porous media for the description of dynamic
hydraulic fracturing, part ii: The crack-opening indicator. Computer Methods in Applied Mechanics and Engineering, 341,
429–442, doi:https://doi.org/10.1016/j.cma.2018.07.006.
URL https://www.sciencedirect.com/science/article/pii/S0045782518303372

[17] Ester, M., H.-P. Kriegel, J. Sander, and X. Xu, 1996: A density-based algorithm for discovering clusters in large spatial
databases with noise. Knowledge Discovery and Data Mining.

[18] Fei, F. and J. Choo, 2020: A phase-field method for modeling cracks with frictional contact. International Journal for
Numerical Methods in Engineering, 121, no. 4, 740–762, doi:https://doi.org/10.1002/nme.6242.
URL https://onlinelibrary.wiley.com/doi/abs/10.1002/nme.6242

[19] Gavagnin, C., L. Sanavia, and L. De Lorenzis, 2020: Stabilized mixed formulation for phase-field computation of deviatoric
fracture in elastic and poroelastic materials. Computational Mechanics, 65, doi:10.1007/s00466-020-01829-x.
URL https://link.springer.com/article/10.1007/s00466-020-01829-x

[20] Geelen, R. J. M., Y. Liu, J. E. Dolbow, and A. Rodríguez-Ferran, 2018: An optimization-based phase-field method for
continuous-discontinuous crack propagation. International Journal for Numerical Methods in Engineering, 116, no. 1, 1–
20, doi:https://doi.org/10.1002/nme.5911.
URL https://onlinelibrary.wiley.com/doi/abs/10.1002/nme.5911

[21] Harville, D. A., 1998: Matrix algebra from a statistician’s perspective. Taylor & Francis.

[22] Heider, Y., 2021: A review on phase-field modeling of hydraulic fracturing. Engineering Fracture Mechanics, 253, 107881,
doi:https://doi.org/10.1016/j.engfracmech.2021.107881.
URL https://www.sciencedirect.com/science/article/pii/S0013794421003143

[23] Hotelling, H., 1933: Analysis of a complex of statistical variables into principal components. Journal of Educational
Psychology, 24, no. 6, 417–441, doi:https://doi.org/10.1037/h0071325.

[24] Hyman, J. D., S. Karra, N. Makedonska, C. W. Gable, S. L. Painter, and H. S. Viswanathan, 2015: dfnworks: A dis-
crete fracture network framework for modeling subsurface flow and transport. Computers & Geosciences, 84, 10–19,
doi:https://doi.org/10.1016/j.cageo.2015.08.001.
URL https://www.sciencedirect.com/science/article/pii/S0098300415300261

[25] Le, D. T., J.-J. Marigo, C. Maurini, and S. Vidoli, 2018: Strain-gradient vs damage-gradient regularizations of softening
damage models. Computer Methods in Applied Mechanics and Engineering, 340, 424–450.

[26] Lee, S., M. F. Wheeler, and T. Wick, 2017: Iterative coupling of flow, geomechanics and adaptive phase-field fracture
including level-set crack width approaches. Journal of Computational and Applied Mathematics, 314, 40–60.
24

[27] Miehe, C., M. Hofacker, and F. Welschinger, 2010: A phase field model for rate-independent crack propagation: Robust
algorithmic implementation based on operator splits. Computer Methods in Applied Mechanics and Engineering, 199, no.
45, 2765–2778, doi:https://doi.org/10.1016/j.cma.2010.04.011.
URL https://www.sciencedirect.com/science/article/pii/S0045782510001283

[28] Miehe, C., S. Mauthe, and S. Teichtmeister, 2015: Minimization principles for the coupled problem of darcy–biot-type
fluid transport in porous media linked to phase field modeling of fracture. Journal of the Mechanics and Physics of Solids,
82, 186–217.

[29] Miehe, C., L.-M. Schaenzel, and H. Ulmer, 2015: Phase field modeling of fracture in multi-physics problems. part i. balance
of crack surface and failure criteria for brittle crack propagation in thermo-elastic solids. Computer Methods in Applied
Mechanics and Engineering, 294, 449–485, doi:https://doi.org/10.1016/j.cma.2014.11.016.
URL https://www.sciencedirect.com/science/article/pii/S0045782514004423

[30] Nguyen, T., J. Yvonnet, Q.-Z. Zhu, M. Bornert, and C. Chateau, 2016: A phase-field method for computational model-
ing of interfacial damage interacting with crack propagation in realistic microstructures obtained by microtomography.
Computer Methods in Applied Mechanics and Engineering, 312, 567–595, doi:https://doi.org/10.1016/j.cma.2015.10.007,
phase Field Approaches to Fracture.
URL https://www.sciencedirect.com/science/article/pii/S0045782515003266

[31] Nicholson, J., 2022: regularizeNd(https://de.mathworks.com/matlabcentral/fileexchange/61436-regularizend). MATLAB

Central File Exchange.
URL https://de.mathworks.com/matlabcentral/fileexchange/61436-regularizend

[32] Osher, S. and R. P. Fedkiw, 2001: Level set methods: an overview and some recent results. Journal of Computational
physics, 169, no. 2, 463–502.

[33] Pauly, M., M. Gross, and L. Kobbelt, 2002: Efficient simplification of point-sampled surfaces. IEEE Visualization, 2002.
VIS 2002., 163–170.

[34] Peerlings, R. d., R. d. Borst, W. d. Brekelmans, J. d. Vree, and I. Spee, 1996: Some observations on localization in non-local
and gradient damage models. European Journal of Mechanics. A, Solids, 15, no. 6, 937–953.

[35] Pijaudier-Cabot, G. and Z. P. Bažant, 1987: Nonlocal damage theory. Journal of engineering mechanics, 113, no. 10,
1512–1533.

[36] Sander, J., M. Ester, H.-P. Kriegel, and X. Xu, 1998: Density-based clustering in spatial databases: The algorithm gdbscan
and its applications. Data Mining and Knowledge Discovery, 2, 169–194.

[37] Santillán, D., R. Juanes, and L. Cueto-Felgueroso, 2017: Phase field model of fluid-driven fracture in elastic media:
Immersed-fracture formulation and validation with analytical solutions. Journal of Geophysical Research: Solid Earth,
122, no. 4, 2565–2589.

[38] Schubert, E., J. Sander, M. Ester, H. P. Kriegel, and X. Xu, 2017: Dbscan revisited, revisited: Why and how you should
(still) use dbscan. ACM Trans. Database Syst., 42, no. 3, doi:10.1145/3068335.
URL https://doi.org/10.1145/3068335

[39] Svendsen, B., P. Shanthraj, and D. Raabe, 2018: Finite-deformation phase-field chemomechanics for
multiphase, multicomponent solids. Journal of the Mechanics and Physics of Solids, 112, 619–636,
doi:https://doi.org/10.1016/j.jmps.2017.10.005.
URL https://www.sciencedirect.com/science/article/pii/S0022509617305495

[40] Tamayo-Mas, E. and A. Rodríguez-Ferran, 2015: A medial-axis-based model for propagating cracks in a regularised bulk.
International Journal for Numerical Methods in Engineering, 101, no. 7, 489–520, doi:https://doi.org/10.1002/nme.4757.
URL https://onlinelibrary.wiley.com/doi/abs/10.1002/nme.4757
 25

[41] Vermilye, J. M. and C. H. Scholz, 1998: The process zone: A microstructural view of fault growth. Journal of Geophysical
Research: Solid Earth, 103, no. B6, 12223–12237, doi:https://doi.org/10.1029/98JB00957.
URL https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/98JB00957

[42] von Wahl, H. and T. Wick, 2022: A high-precision framework for phase-field fracture interface reconstructions with
application to stokes fluid-filled fracture surrounded by an elastic medium. arXiv preprint arXiv:2212.07982.

[43] Wu, J.-Y. and W.-X. Chen, 2022: On the phase-field modeling of fully coupled chemo-mechanical deterioration
and fracture in calcium leached cementitious solids. International Journal of Solids and Structures, 238, 111380,
doi:https://doi.org/10.1016/j.ijsolstr.2021.111380.
URL https://www.sciencedirect.com/science/article/pii/S0020768321004467

[44] Yang, L., Y. Yang, H. Zheng, and Z. Wu, 2021: An explicit representation of cracks in the variational phase
field method for brittle fractures. Computer Methods in Applied Mechanics and Engineering, 387, 114127,
doi:https://doi.org/10.1016/j.cma.2021.114127.
URL https://www.sciencedirect.com/science/article/pii/S0045782521004588

[45] Yoshioka, K., D. Naumov, and O. Kolditz, 2020: On crack opening computation in variational phase-field models for
fracture. Computer Methods in Applied Mechanics and Engineering, 369, 113210.

[46] Zeng, J., M. Zhang, E. Yang, F. Tian, and L. Li, 2022: A tracking strategy for multi-branched crack tips in phase-
field modeling of dynamic fractures. International Journal for Numerical Methods in Engineering, 123, no. 3, 844–865,
doi:https://doi.org/10.1002/nme.6879.
URL https://onlinelibrary.wiley.com/doi/abs/10.1002/nme.6879

[47] Zhang, J., H. Yu, W. Xu, C. Lv, M. Micheal, F. Shi, and H. Wu, 2022: A hybrid numerical approach for hydraulic fractur-
ing in a naturally fractured formation combining the xfem and phase-field model. Engineering Fracture Mechanics, 271,
108621.

[48] Ziaei-Rad, V., L. Shen, J. Jiang, and Y. Shen, 2016: Identifying the crack path for the phase field approach to
fracture with non-maximum suppression. Computer Methods in Applied Mechanics and Engineering, 312, 304–321,
doi:https://doi.org/10.1016/j.cma.2016.08.025, phase Field Approaches to Fracture.
URL https://www.sciencedirect.com/science/article/pii/S0045782516310209

View publication stats