Efficient Crack Length Measurement Using A* Shortest Path Methodology for a Phase-Field Fracture Framework

Abstract

Accurately measuring a crack length is a crucial aspect of experimental fracture tests. In this work, we present an innovative application of the A* (A-star) shortest path methodology to track different shapes of cracks from numerical simulations. This approach is highly efficient, significantly improving the speed and accuracy of crack length measurements. Furthermore, we introduce a modified weight cost function that follows the crack path in the damage field, enhancing the accuracy of our method. The effectiveness of the proposed procedure is shown by fabricating damage fields with different geometry and good agreement when compared to the exact values. In addition, we evaluate a time-dependent crack propagation case, achieving high accuracy. We present all features and steps of the procedure to showcase its efficacy in accurately measuring the length of a crack path. Finally, we validate our method using a phase-field fracture framework and compare it with the compliance technique. The results show that the proposed method is applicable in finite element analyses with recovering accurate results.

Keywords
Phase Field; Crack length; A-Star; Shortest-Path

Graphical Abstract

1 INTRODUCTION

The fracture of materials has been of concern in the industry and academia over the years. In response to this issue, experimental studies have been conducted to better understand this phenomenon and develop strategies to mitigate catastrophic failures.

The introduction and subsequent widespread utilization of material parameters and their correlation, such as fracture toughness, $K$, and crack propagation rate, $Da/DN$ (obtained from the Paris curve), have made the investigation of crack behavior simpler and the prediction of potential mechanical failure of materials. These parameters are closely linked to crack length, for which precise measurement is crucial.

Numerous crack length measurement methods are available for data obtained from experimental tests. Examples include the compliance technique (Saxena and Hudak, 1978; Kanters et al., 2015), the direct length measurement using an optical microscope (Lee et al., 2015), and the direct current potential drop (Bär, 2020). These methods rely on force and displacement measurements, which require accurate experimental equipment. However, alternative methods have been developed to address this challenge that requires no experimental output. These methods include thermoelasticity (Díaz et al., 2004), photoelasticity (Jobin et al., 2020), and digital image correlation (DIC) (Mokhtarishirazabad et al., 2018).

Researchers have developed new methods and improvements for measuring crack length over the years. In a study conducted by Creel et al. (2009), compliance parameters for fracture testing of an anisotropic biological material were calibrated. Fünfschilling et al. (2010) compared crack length measurements obtained using optical and compliance methods when applied to R-curves. In more recent years, Yuan et al. (2021) utilized artificial intelligence based on convolution neural networks (CNN) and digital image processing to monitor crack length. Similarly, Farahani et al. (2022) combined electronic speckle pattern interferometry (ESPI) with digital image correlation (DIC) to obtain the displacement field, calculate the stress intensity factor (SIF), and monitor the crack opening evolution.

Recently, phase field-based methods have been developed for modeling fracture phenomena. They use a continuous damage variable, usually denoted by $\phi$ with range from $0$ for undamaged material and $1$ for fully fractured material, to represent physical crack discontinuity. Many publications have reported good qualitative and quantitative agreement with experimental tests. However, measuring crack length in numerical simulations using a predetermined threshold damage remains challenging. Klinsmann et al. (2015), Zhang et al. (2018), Carrara et al. (2020), and Liu et al. (2020) used a damage threshold of $\mathrm{\phi }=0.95$ and emphasize that the uncertainty in measuring the crack length contributes significantly to the discrepancy of results. This method becomes inaccurate when dealing with general shape cracks under mixed-mode loading (Zhang et al., 2017) or dynamic crack propagation cases with crack branching (Bleyer et al., 2017). For accurate numerical analysis, better numerical methodologies are required to follow experimental tests and evaluate crack propagation (ASTM E647, 2016) and fracture toughness (ASTM E399, 2017). Recently, Gehri et al. (2020) introduced innovative techniques that use digital image correlation (DIC) and image processing to measure cracks experimentally, resulting in improvements in multiple crack detection. Nonetheless, it is important to note that the use of image filters can substantially increase computational costs, which is dependent on the analysis periodicity and the total number of image pixels.

A-Star ($A^{*}$) is a popular shortest-path method with numerous applications, including map routing (Wang et al., 2014), autonomous vehicle path planning (Erke et al, 2020), movement of military troops (Dawid and Pokonieczny, 2021), and commercial video games (Lawande et al, 2022). This method was introduced by Peter Hart, Nils Nilsson, and Bertram Raphael (Hart et al, 1968) and can be seen as an extension of Dijkstra's algorithm (Dijkstra, 1959). Considering heuristics in$A^{*}$ leads to improved performance and efficiency in finding the shortest path with low cost and fast implementation.

This work proposes a novel application of the $A^{*}$ pathfinding algorithm to numerically measure the crack length using the finite element method (FEM) and a phase field fracture model. The key concept is to adapt the $A^{*}$ method, typically associated with graph analysis, to the finite element mesh and use the damage scalar field as a guide to track the crack path. For this purpose, the vertices, edges, and faces of graphs are associated, respectively, to the nodes, edges, and elements of the finite element meshes to make possible the adequate operation of the $A^{*}$ method.

This paper is structured as follows. Firstly, we provide a concise introduction to the phase field fracture model and the equations system to be solved. Subsequently, we present the typical $A^{*}$ pathfinding method, followed by the novel weighted cost function and the procedure for measuring the crack length. Finally, we illustrate the effectiveness of our approach by presenting the path length results for artificially fabricated damage fields and a crack resulting from fracture simulation with the phase field methodology.

2 PHASE FIELD MODELLING

Phase field methodology utilizes the continuous variable $\phi$ (see Figure 1) to describe an interface problem. It was initially introduced to solidification (Boettinger et al., 2002) and material phase transformation (Penrose and Fife, 1990). More recently, this methodology has found applications in material modelling for fracture analysis, fatigue (Boldrini et al, 2016; Alessi et al, 2018 and Haveroth et al, 2020), anisotropy (Petrini et al, 2021), and hydrogen embrittlement (Martínez-Pañeda et al, 2018).

Figure 1
Damage phase field where $\mathrm{\phi }\approx 1$ represents a physical crack.

A phase field model with thermodynamic consistency is derived using the balance equations of continuum mechanics, resulting in a general and robust framework. In this work, the formulation is given according to Boldrini et al. (2016) and Haveroth et al. (2020).

The conservation of mass is given by

$\dot{\rho }+\rho div\left(\mathit{v}\right)=0$,(1)

where $\mathrm{\rho }$ is the volumetric mass density and $\mathit{v}$ is velocity vector field given by time derivative of the displacement field, $\mathit{u}$, both represented in the Eulerian description.

The principle of virtual power (PVP) is applied to obtain the following equilibrium equation at the macro scale for an arbitrary part ${\mathcal{D}}_{\mathcal{t}}$ for the body $\mathcal{B}$ at time $t$:

$\left\{\begin{array}{c}\rho \dot{\mathit{v}}=div\left(\mathit{T}\right)+\rho fin{\mathcal{D}}_{\mathcal{t}}\\ Tn=tin\partial {\mathcal{D}}_{\mathcal{t}}\end{array}\right.$,(2)

where $\mathit{T}$ is the Cauchy stress tensor field, $\mathit{f}$ is the body load vector field and $\mathit{n}$ is the normal vector field to the boundary $\partial {\mathcal{D}}_{\mathcal{t}}$.

The PVP is also used at the microscale resulting in the following equilibrium equation:

$\left\{\begin{array}{c}div\left(\mathit{h}\right)-b+\rho a=0in{\mathcal{D}}_{\mathcal{t}}\\ h\cdot n=t_{h\mathrm{}}in\partial {\mathcal{D}}_{\mathcal{t}}\end{array}\right.$,(3)

where $h$ is the microstress vector field, $b$ is the microforce scalar field, and $a$ is the density of external body microforce.

The balance of energy, or the first principle of thermodynamics, is given by

$\rho \dot{e}=\mathit{T}:D+b\dot{\phi }+\mathit{h}\cdot \nabla \dot{\phi }-div\left(\mathit{q}\right)+\rho rin{\mathcal{D}}_{\mathcal{t}}$,(4)

where $e$ is the specific internal energy density, $\mathit{q}$ is the heat flux vector field, and $r$ is the specific heat source density.

Finally, the second law of thermodynamics, stated by the generalized Clausius-Duhem inequality, is

$\rho \dot{\eta }\ge -div\left(\frac{\mathit{q}}{\theta }+\mathit{k}\right)+\rho \left(\frac{r}{\theta }+\omega \right)in{\mathcal{D}}_{\mathcal{t}}$,(5)

where $\mathrm{\eta }$ is the specific entropy density, $\mathit{k}$ is the entropy flux due to microscopic evolution, $\mathrm{\theta }$ is the absolute temperature, $r/\mathrm{\theta }$ is the specific entropy production, and $\mathrm{\omega }$ is the entropy production due to microscopic evolution.

To obtain thermodynamically consistent relations, the constitutive properties of the material are expressed in terms of the free-energy potential, $\mathrm{\psi }$, and the pseudo-potential of dissipation, ${\mathrm{\psi }}_d^n$. Using the balance of energy of Eq. 4, the entropy inequality, Eq. 5, is rewritten in terms of the specific free-energy $\mathrm{\psi }=e-\mathrm{\theta }\mathrm{\eta }$ by

$-\rho \left(\dot{\psi }+\eta \dot{\theta }\right)+\mathit{T}:\mathit{D}+b\dot{\phi }+\mathit{h}\cdot \nabla \dot{\phi }-\frac{1}{\theta }\mathit{q}\cdot \nabla \theta +\theta div\left(\mathit{k}\right)-\rho \omega \theta \ge 0$.(6)

The free energy potential depends on the internal variables as

$\left\{\begin{array}{c}\psi =\psi \left(\mathrm{\Gamma }\right)\\ \Gamma =\left(\rho ,\theta ,\phi ,\nabla \rho ,\nabla \theta ,\nabla \phi ,\mathit{E}\right)\end{array}\right.$.(7)

Moreover, the terms $\mathit{T}$, $b$, $\mathit{h}$ and $\mathit{q}$ are split in their reversible (non-dissipative) and irreversible (dissipative) parts and represented as ${\left(\cdot \right)}^r$and ${\left(\cdot \right)}^{ir}$ respectively. In addition, it is assumed that ${\mathit{h}}^{\left(ir\right)}\equiv 0$ and ${\mathit{q}}^{\left(r\right)}\equiv 0$.

For the sake of brevity, some manipulations from the previous equations will be omitted here. Further details can be found in Boldrini et al. (2016).

Reversible terms are chosen such that, for any admissible process, there is no entropy increase. Considering irreversible terms, and after some simplifications, we have the following equations:

$b^{\left(ir\right)}\dot{\phi }+{\mathit{T}}^{\left(ir\right)}:\mathit{D}-\frac{1}{\theta }{\mathit{q}}^{\left(ir\right)}\ge 0$.(8)

To satisfy the previous inequality, we define the pseudo-potential of dissipation by

${\psi }_d^n={\psi }_d^n\left(\dot{\phi },\mathit{D},\nabla \theta ,\tilde{\mathrm{\Gamma }}\right)$ (9)

with

$\tilde{\mathrm{\Gamma }}=\left(\rho ,\theta ,\phi ,\nabla \rho ,\mathit{E}\right)$.(10)

If ${\psi }_d^n\left(\dot{\mathrm{\phi }},\mathit{D},\nabla \mathrm{\theta },\tilde{\mathrm{\Gamma }}\right)\ge 0,\forall \left\{\dot{\mathrm{\phi }},\mathit{D},\nabla \mathrm{\theta },\tilde{\mathrm{\Gamma }}\right\}$, ${\mathrm{\psi }}_d^n\left(\mathrm{0,0},0,\tilde{\mathrm{\Gamma }}\right)\ge 0$ and ${\mathrm{\psi }}_d^n$ is continuous and convex with respect to the variables $\dot{\mathrm{\phi }},\mathit{D}$ and $\nabla \mathrm{\theta }$, Eq. 8 is satisfied.

Collecting the previous equations and after some simplifications, the system of governing equations is expressed by

$\left\{\begin{array}{c}\dot{\rho }+\rho div\left(\mathit{v}\right)=0\\ \dot{\mathit{u}}=v\\ \rho \dot{\mathit{v}}=div\left(\mathit{T}\right)+\rho f\\ T=\rho {\partial }_{\mathit{E}}\psi -{\rho }^2{\partial }_{\rho }\psi I-\rho sym\left(\nabla \mathcal{F}\otimes {\partial }_{\nabla \mathcal{F}}\psi +\nabla \phi \otimes {\partial }_{\theta }\psi \right)+{\partial }_{\mathit{D}}{\psi }_d^n\\ {\partial }_{\dot{\phi }}{\psi }_d^n=div\left(\theta {\partial }_{\nabla \theta }{\psi }_d^n\right)+T:D+\left(\rho {\partial }_{\phi }\psi +{\partial }_{\dot{\phi }}{\psi }_d^n\right)\dot{\phi }+\rho {\partial }_{\nabla \phi }\psi \cdot \nabla \dot{\phi }+\rho r\\ e=\psi -\theta {\partial }_{\theta }\psi \end{array}\right.$.(11)

This system of equations is very general. Changing the free energy $\mathrm{\psi }$ and the pseudo-potential ${\mathrm{\psi }}_d^n$, different thermodynamically consistent models that account for fracture phenomena can be obtained for many materials.

Let us consider the particular case of an isothermal, isotropic, and linear elastic material with damage.

The volumetric free-energy density for a material subjected to damage is given by the sum of the elastic energy density $\mathcal{E}$ and the energy density related to damage $\mathcal{J}$ as

${\rho }_0\psi \left(\phi ,\nabla \phi ,\mathit{E}\right)=\mathcal{E}+\mathcal{J}$,(12)

where a nearly-incompressible material is such that

$\rho \left(\mathit{x},t\right)={\rho }_0$.(13)

The elastic strain energy density of the damaged material is obtained from the respective energy of the undamaged state by

$\mathcal{E}=g_e{\mathcal{E}}_0$.(14)

where $g_e$ is the degradation function and ${\mathcal{E}}_0$ is the elastic strain energy density of the undamaged material. The degradation function may depend on several variables. In this work, we assume that

$g_e\left(\phi \right)={\left(1-\phi \right)}^{\kappa }$,(15)

where $\mathrm{\kappa }$ is a material dependent real constant.

Moreover, we split the elastic strain energy density, as proposed by Amor et al. (2009), into its volumetric and deviatoric (volume-preserving) contributions such that

$\mathcal{E}\left(\mathit{E},\phi \right)=g_e^{+}{\mathcal{E}}_0^{+}\left(\mathit{E}\right)+g_e^{-}{\mathcal{E}}_0^{-}\left(\mathit{E}\right)$,(16)

where

${\mathcal{E}}_0^{+}=\frac{1}{2}K{⟨tr\left(\mathit{E}\right)⟩}_{+}^2+G\left({\mathit{E}}_{\text{dev}}:{\mathit{E}}_{\text{dev}}\right)$,(17)

and

${\mathcal{E}}_0^{-}=\frac{1}{2}K{⟨tr\left(\mathit{E}\right)⟩}_{-}^2$,(18)

where $G$ and $K$ are the shear and the bulk modulus of the material, respectively. The operator ${⟨\cdot ⟩}_{\pm }$ is defined as ${⟨\mathrm{\xi }⟩}_{\pm }=\left(\mathrm{\xi }\pm \left|\mathrm{\xi }\right|\right)/2$.

The crack energy density $\mathcal{J}$ is the sum of the energy density associated with the transition layers and the total bulk energy density. Therefore,

$\mathcal{J}\left(\phi ,\nabla \phi \right)=G_c\left(\frac{\gamma }{2}{\left|\nabla \phi \right|}^2+\frac{{\phi }^2}{2\gamma }\right)$,(19)

where $\mathrm{\gamma }$ is a parameter related to the width of the fracture layers and $G_c$ is the Griffith fracture energy parameter.

The non-plastic pseudo-potential of dissipation is defined by

${\psi }_d^n\left(\dot{\phi },\nabla \theta ,\tilde{\mathrm{\Gamma }}\right)=\frac{1}{2}\tilde{\lambda }\left(\tilde{\mathrm{\Gamma }}\right){\left|\dot{\phi }\right|}^2$,(20)

where $\tilde{\mathrm{\lambda }}$ controls the damage change rate, see Boldrini et al. (2016).

A quasi-static process is here considered for the motion equation. General thermodynamically consistent equations are expressed in the following way:

$\left\{\begin{array}{c}0=\frac{1}{{\rho }_0}div\sigma -\frac{\gamma G_c}{{\rho }_0}div\left(\nabla \mathrm{}\phi \mathrm{}\otimes \mathrm{}\nabla \mathrm{}\phi \right)+f\\ \dot{\phi }=\frac{\gamma G_c}{\tilde{\lambda }}\Delta \phi -\frac{1}{\tilde{\lambda }}\left({\partial }_{\phi }{g}_e^{+}{\mathcal{E}}_0^{+}+{\partial }_{\phi }{g}_e^{-}{\mathcal{E}}_0^{-}\right)-\frac{1}{\tilde{\lambda }\gamma }{G}_c\phi \end{array}\right.$,(21)

with

$\sigma ={\partial }_{\mathit{E}}\mathcal{E}$,(22)

$\frac{1}{\tilde{\lambda }}=\frac{\tilde{c}}{\left(1-\phi +\delta \right)},$ (23)

where $\tilde{c}$ is a constant that controls the damage growth rate and $\mathrm{\delta }$ a small constant to avoid numerical singularity.

3 DETERMINATION OF CRACK PATH

This section describes the proposed procedure for tracking crack paths using a modified version of the $A^{*}$ method.

The $A^{*}$ based methods require prior knowledge of the starting and ending points of the path. While the crack nucleation position $\mathcal{A}$ is generally known from experimental tests, finding the respective initial node of the mesh can be challenging in complex geometries. However, it can be estimated using the damage phase field variable and the maximum von Mises stress. Once $\mathcal{A}$ is found, it remains constant along the procedure to determine the crack path.

Similarly, the crack tip $\mathcal{B}$ is chosen to be located at the point with the maximum local von Mises stress. Its position must be updated at every load step of the simulation to track the crack path.

3.1 Weighted $A^{*}$ Method

The $A^{*}$ algorithm is a commonly used informed search approach in pathfinding problems. The algorithm aims to find the path from a starting point $\mathcal{A}$ (crack nucleation site) to a goal point $\mathcal{B}$ (crack tip) with the minimum value of the cost function $f$ for moving between $\mathcal{A}$ and $\mathcal{B}$, which are both known positions.

The path $\mathcal{P}$ is defined as a set of consecutive nodes, denoted by $n$, in the finite element mesh, where $n$ is a positive integer. The path $\mathcal{P}$ connects $\mathcal{A}$ and $\mathcal{B}$ with the minimum cost by selecting the nodes with the least accumulated cost.

The shortest path between points $\mathcal{A}$ and $\mathcal{B}$ are described in terms of mesh nodes denoted by the set $S_{\mathcal{A}{n}^p}$ of parent nodes which start at point $\mathcal{A}$. Variables with superscript ${\left(\cdot \right)}^n$ denote values computed at node $n$, while those variables with superscript ${\left(\cdot \right)}^p$ represent the corresponding values for a parent node.

A cost function is typically used at each node $n$ and expressed as

$f^n=g^n+h^n$.(24)

Term $g^n$ represents the accumulated distances between the previous parent nodes and the distance between node $n$ and its parent node $n^p$; term $h^n$ is an estimated heuristic cost that approximates the distance between node $n$ and the final goal point $\mathcal{B}$ (Tortora and Francese, 2003).

Specifically, $g\left(n\right)$ can be expressed as $g\left(n\right)=S_{\mathcal{A}{n}^p}+{\left|\left|n^p-n\right|\right|}_2$, where $$\left|\right|\cdot |{|}_2$ is the Euclidean norm. Similarly, $h\left(n\right)={\left|\left|\mathcal{B}-n\right|\right|}_2$.

Data for each considered node is allocated in the structure called OPEN, which stores the node number and its parent node number. The cost function $g$ is recursively calculated using Algorithm 1, while the cost function $h$ is calculated using the Euclidean distance mentioned before.

Algorithm 1: Recursive function for $\mathit{g}$

Initialize: $g=0$
Input: $OPEN$, $g$, $n$, $Mesh$
Output: $g,n$
1:	Identify the parent node of the actual node in $OPEN$
2:	Calculate: $g=g+{\left|\left|n-n^p\right|\right|}_2$
3:	if Current node is the start node, then
4:	Return to the main program
5:	end if
6:	Update: $n=n^p$
7:	Call this function

The method involves tracking nodes with the minimum cost from point $\mathcal{A}$ to $\mathcal{B}$ in order to find the shortest path between them. At each iteration, the program allocates information for the neighboring nodes $n_1$, $n_2$, $\dots$, $n_{\mathrm{\eta }}$ for the current node $n$ also in the OPEN structure. These neighboring nodes can be chosen from the adjacent nodes (i.e., nodes connected by an element edge) or the nodes that belong to the neighboring elements (see Figure 2). When searching for adjacent nodes in a quadrilateral element, the path length will depend on both the element edges and the crack directions.

Figure 2
Finite element mesh with a central node $n$ (black) and its neighboring nodes connected by an edge (red) and the nodes of the neighboring elements (blue).

Node $n$ is appended to the CLOSED matrix, which contains the nodes already analyzed. The OPEN matrix must be sorted in ascending order of the cost function $f$ before the new iteration. The loop criterion is achieved when the analyzed node is the goal node. The $A^{*}$ procedure is summarized in the Algorithm 2.

The nodes in CLOSED matrix are sequenced backward for building the nodal path sequence. For example, starting from the goal node (the first position) then the parent node provides the next node in the series, and so forth, up to the start node.

Algorithm 2: Weighted ${\mathit{A}}^{\mathit{*}}$

Input: Mesh, damage field
Output: Nodal path
1:	Allocate the start node at $OPEN$
2:	while Node $\ne$ Goal Node, do
3:	for Every neighbor node by edge or by element, do
4:	if node is not in $OPEN$ or not in $CLOSED$, then
5:	Calculate $g$, $h$, and $f$
6:	Store the node information in $OPEN$
7:	end if
8:	end for
9:	Store Node in $CLOSE$
10:	Remove Node From $OPEN$
11:	Sort $OPEN$ in ascending order of cost function $f$
12:	end while

The previous procedure is efficient in searching for the shortest path in a finite element mesh. However, it does not guarantee that the path will pass through the crack. To enforce this condition, the cost functions $g$ and $h$ are weighted by the damage parameter, $\mathrm{\phi }\left(n\right)$, as

$f\left(n\right)=w\left(\phi \left(n\right)\right)\left(g\left(n\right)+h\left(n\right)\right)$.(25)

where $w\left(\mathrm{\phi }\left(n\right)\right)$ is the weight function, which reduces the cost functions where the nodal damage, ${\mathrm{\phi }}^n$, approaches $1$, as indicated by the considered degradation function

$w\left(\phi \left(n\right)\right)={\left(1-\phi \right)}^{\alpha }$ (26)

where $\mathrm{\alpha }\in {\mathbb{R}}^{+}$.

There are different weight functions, $w\left(\mathrm{\phi }\left(n\right)\right)$, proposed in the literature. In this case, it is crucial to ensure that the total cost decreases as the damage increases, forcing the path to pass through the most damaged nodes. Consequently, the weighted cost function no longer searches for the globally shortest path. It should be noted that the classical $A^{*}$ method is recovered when $\mathrm{\phi }=0$ or $\mathrm{\alpha }=0$.

3.2 Approximation of the crack length

We propose a simple and effective method for calculating the length of the crack path based on a line integral.

The parametrization of the curve $C$ in terms of the local coordinate $\mathrm{\xi }\in \left[-\mathrm{1,1}\right]$ from one node with coordinates $\left[x_i,y_i\right]$ to another node with coordinates $\left[x_f,y_f\right]$ is given by

$I={\int }_{-1}^{1}f\left(x\left(\xi \right),y\left(\xi \right)\right)\left|{\mathit{r}}^{\mathrm{\text{'}}}\left(x\left(\xi \right),y\left(\xi \right)\right)\right|d\xi$ (27)

where

$\left\{\begin{array}{c}{\mathit{r}}^T\left(\xi \right)=\left[x\left(\xi \right)y\left(\xi \right)\right]\\ x\left(\xi \right)=\frac{x_i+x_f}{2}+\frac{x_f-x_i}{2}\\ y\left(\xi \right)=\frac{y_i+y_f}{2}+\frac{y_f-y_i}{2}\end{array}\right.$ (28)

In the finite element mesh, the scalar field $\mathrm{\phi }$ is interpolated using the one-dimensional linear nodal shape functions $N_i$ expressed in the local coordinates. Eq. 27 is then rewritten for an element that contains two consecutive nodes of the path $\mathcal{P}$ as

$I=\int N_i\left(x\left(\xi \right),y\left(\xi \right)\right){\phi }_i\left|{\mathit{r}}_j^{\mathrm{\text{'}}}\left(x\left(\xi \right),y\left(\xi \right)\right)\right|d\xi$ (29)

Since ${\mathit{r}}_j\mathrm{\text{'}}\left(x\left(\mathrm{\xi }\right),y\left(\mathrm{\xi }\right)\right)$ is independent of $\mathrm{\xi }$, Eq.29 is rewritten as

$I\mathrm{}=\mathrm{}\underset{\text{L}}{\underbrace{\left|\mathit{r}{\mathrm{\text{'}}}_j\left(x\left(\xi \right),y\left(\xi \right)\right)\right|}}\mathrm{}\underset{\text{M}}{\underbrace{\int N_i\left(x\left(\xi \right),y\left(\xi \right)\right){\phi }_{i}d\xi }}$ (30)

where $L$ is the curve length and $M$ is the nodal damage weighted by the nodal shape function.

Finally, the crack length is given by

$a_{A^{\mathrm{*}}}=\sum _{}^{\mathcal{P}}L\cdot M$ (31)

For elements with $M\ge {\mathrm{\phi }}_t$, i.e., almost fully cracked, $$M$ is considered $1$.

4 RESULTS

This section presents fabricated solution examples to illustrate the capabilities and features of the proposed method. Additionally, a time-dependent case that utilizes the phase-field crack solution is presented to validate the effectiveness of the method.

4.1 Fabricated damage fields

To illustrate the influence of the weight function and validate the proposed $A^{*}$ method, the following fabricated damage field is considered:

$\phi \left(\mathit{x},0\right)=max\left(\phi \left(\mathit{x},0\right),0.99\cdot \mathrm{e}\mathrm{x}\mathrm{p}\{-1\cdot {10}^2{\left(x-x_{\phi }\right)}^2\}\cdot \mathrm{e}\mathrm{x}\mathrm{p}\{-1\cdot {10}^2{\left(y-y_{\phi }\right)}^2\}\right)$ (32)

where $x_{\mathrm{\phi }}$ and $y_{\mathrm{\phi }}$ are the coordinates of the nodes at a critical damage state.

We analyze the two crack configurations of Figure 3 to test the effectiveness of the proposed method. The first configuration, $C_1$ (Figure 3a), represents a curved crack configuration, while the second configuration $C_2$ (Figure 3b), depicts a branched crack configuration. Circular shapes were chosen for the crack configurations motivated by the fact that approximating circular geometry using straight lines is particularly challenging. Consequently, obtaining good results for these cases ensures that simpler cases can also be accurately approximated.

Figure 3
Two-dimensional domains: (a) $C_1$ and (b) $C_2$, illustrating the starting points (blue) and goal points (red) associated with each domain.

The critical damage nodes coordinates are defined as

$C_1:{\left(r-\delta \right)}^2\le \left(x_{\phi }^2+y_{\phi }^2\right)\le {\left(r+\delta \right)}^2$ (33)

$C_2:\left\{\begin{array}{c}{x}_i\le x_{\phi }\le x_f\\ y_i\le y_{\phi }\le y_f\\ {\left(r-\delta \right)}^2\le \left(x_{\phi }^2+y_{\phi }^2\right)\le {\left(r+\delta \right)}^2\end{array}\right.$ (34)

where $r$ is the radius of the circumference and $\mathrm{\delta }$ is a small increment.

The starting and goal points coordinates ($x_{\mathrm{\phi }},y_{\mathrm{\phi }}$) are, respectively, $\left(\mathrm{0.0,0.6}\right)$ and $\left(\mathrm{0.6,0.0}\right)$ for $C_1$. The geometry $C_2$ represents a straight crack followed by branching and described by two paths. The starting point coordinate is $\left(\mathrm{0.0,0.5}\right)$ and the two goal points coordinates are $\left(\mathrm{1.0,0.8}\right)$ and $\left(\mathrm{1.0,0.2}\right)$.The mesh domain is defined in the interval $\left[01\right]$ for both $x$ and $y$ coordinates and discretized using $m$ equally spaced elements in each direction.

The resulting paths for two different choices of $w\left(\mathrm{\phi }\right)$ are presented in Figure 4 and 5 using Algorithms 1 and 2. Firstly, the usual $A^{*}$ method is recovered for $w\left(\mathrm{\phi }\right)=1$ and the shortest path ${\mathcal{P}}_{\mathcal{s}}$ using Eq. 24 is obtained and illustrated in Figs. 4b and 5b. For the weighted $A^{*}$ method using Eq. 26, the path ${\mathcal{P}}_{\mathcal{w}}$ pass through the nodes with the largest values of $\mathrm{\phi }$ as shown in Figs. 4c and 5c.

Figure 4
(a) Scalar field for the $C_1$ curve and the obtained paths associated to the weigh function $g$ where (b) $\alpha =0$ and (c) $\alpha =1$.

Figure 5
(a) Scalar field for the $C_2$ curve and the obtained paths associated to the weigh function $g$ where (b) $\alpha =0$ and (c) $\alpha =1$.

The proposed method showed qualitative agreement compared to the original shape. However, a quantitative analysis must be considered to validate the methodology. For this purpose, the method presented in Section 3.2 is applied, and the crack length is obtained for different mesh refinements ($m=\mathrm{50,100,200}$) and different neighbor search options (by edge or by element).

The results are summarized in Table 1 and the accuracy depends on the path shape and on the choice of searching for neighbor nodes (by edge or by element). Evidently, a circular shape in a mesh with squared elements will present lower accuracy when compared to a straight line. In this case, searching for neighbors by node provides more possible nodes and directions improving the path length result.

Table 1
Convergence analysis for the examples with fabricated damage fields.

Geometry	Exact Length	Path Length [m] (% error)
		$m=50$		$m=100$		$m=200$
		By edge	By element	By edge	By element	By edge	By element
${\mathit{C}}_1$	0.942478	0.80	0.989	1.200000	0.989117	1.200000	0.986188
		(15.12%)	(4.95%)	(27.32%)	(4.64%)	(27.32%)	(4.64%)
${\mathit{C}}_2$	1.171239	1.294416	1.182843	1.300000	1.188701	1.300000	1.185772
		(10.52%)	(0.99%)	(10.99%)	(1.49%)	(10.52%)	(1.24%)

4.2 Crack under propagation

The proposed method has a practical application in tracking the path of a stable crack over time. Specifically, the crack path is measured at regular intervals during the simulation, between the initial time $t_i$ and final time $t_f$. This technique is commonly used in fatigue crack propagation tests to determine the propagation velocity of a stable crack under a given stress intensity factor. The resulting data is represented by the Paris-Erdogan curve, which provides valuable material properties, including the threshold stress intensity factor and stable propagation parameters.

The considered square domain with side $1m$ is discretized into structured meshes with $m^2$ elements. Again, the damage scalar field is given by Eq. 32 and updated every time step. The simulation time ranges from $t_i=0$ to $t_f=50s$ with time step $∆t=1s$.

A straight crack pattern is here considered. The nodal positions are constant for every time step and the critical damaged nodal coordinates are given by the curve

$C_3:\left\{\begin{array}{c}{x}_i\le x_{\phi }\le x_f\left(t\right)\\ y_i\le y_{\phi }\le y_f\end{array}\right.$ (35)

with

$x_f\left(t\right)=x_{fi}+\left(L_x-x_{fi}\right)t^4$ (36)

where $x_f\left(t\right)$ is the time-dependent crack tip position function, $x_{fi}$ is the starting crack tip coordinate in the $x$-axis and $L\_x$ is the width of the domain in $x$-direction. The starting crack parameters are $x_{fi}=0.7m$, $x_i=0m$, $y_i=y_f=0.5m$.

The crack length is measured using three meshes with different refinements ($m=50,100,200$) and compared with the exact value given by Eq. 35 and shown in Figure 6.

Figure 6
Crack length measured over time for different mesh refinements: (a) $m=50$, (b) $m=100$ and (c) $m=200$

The proposed method can closely follow the damage field and the crack tip, even with a coarser mesh. Evidently, the accuracy of the result depends on the mesh refinement. For the coarsest mesh, Figs. 6a-6b shows that the crack length evolution behaves as a step function where the damage rate is lower compared to Fig. 6c. Therefore, in the crack propagation simulation, the result accuracy is dependent on the measurement frequency ($f$), the crack speed ($v$), and the mesh refinement ($r$) in such a way that $f\propto v,r^{-1}$.

The average percentage error is calculated for the three mesh refinements with values of $0.92\%,0.51\%$, and $0.27\%$, respectively.

To avoid dependence on the sampling frequency, the measurements are counted only when the crack advances to a new node. Figure 7 shows the results for $m=50$. The differences are more pronounced in the beginning but more accurate when measured. The average percentage error for this case is $0.62\%$.

Figure 7
Crack length measured over time for $m=50$ excluding repeated measurements.

4.3 Comparison with the compliance method

The resulting crack length obtained from the $A^{*}$ method is compared with the crack length using the compliance technique. The compliance technique is often used in experimental tests since it is simpler to obtain the crack length when compared with other more accurate methods. The length of a straight crack can be indirectly measured using the crack mouth opening displacement (CMOD) and the applied force. For a standard specimen geometry, the crack length results from a calibrated function of the CMOD, the applied force, $P$, and the effective stiffness, $E\text{'}$. Also, the geometric parameters $W$ and $B$ are, respectively, the width (measured from the center of the fixation holes) and the thickness. The C(T) specimen is used with dimensions shown in Figure 8 according to the ASTM E399 standard (ASTM E399, 2017).

Figure 8
Comparison between the compliance and $A^{\mathrm{*}}$ crack measurement methods.

A fracture toughness test (mode I opening test) in the specimen under displacement control is simulated using the phase-field model stated in Eq. 21. The CMOD and the reaction force are provided from the finite element analysis, which are inputs for the crack length functional given by

$a_c=W\left(1.0-4.5\cdot U+13.157\cdot U^2-172.551\cdot U^3+879.944\cdot U^4-1514.671\cdot U^5\right)$ (37)

with

$U={\left(\sqrt{E^{\mathrm{\text{'}}}CB}+1\right)}^{-1}$ (38)

where $a_C$ is the compliance crack length, $U$ is a dimensionless variable, and $C$ is the secant compliance given by

$C=\frac{CMOD}{P}$.(39)

Eq. 37 is valid for crack lengths of $0.2W\le a_c\le 0.8W$ (ASTM E399, 2017).

For the phase-field simulations, the parameters are $\tilde{c}=5.0\times {10}^{-5}$, $E=44GPa$, $\mathrm{\nu }=0.3$, $\mathrm{\Delta }t=0.0167s$, $G_c=1830J/m^2$ and $\mathrm{\gamma }=3.5\times {10}^{-3}m$. The analysis is done under the elastic regime, without plastic deformations. In this study, we have discretized the geometry utilizing over 75000 linear quadrilateral elements which have been strategically concentrated along the expected crack path.

The crack length is computed using the proposed $A^{*}$ method every 20 time steps and also with the compliance technique every time step with $E^{\text{'}}=53GPa$. The mesh is unstructured and refined in the expected crack path with an average element width of $5\times {10}^{-5}m$. The goal node is located at the maximum von Mises stress and a physical crack is summed when $\mathrm{\phi }\ge 0.99$.

The crack length measurements for the proposed $A^{*}$ and the compliance methods are shown in Figure 9. Notably, the weighted $A^{*}$ approach demonstrates remarkable agreement with the compliance method across the entire crack propagation.

Figure 9
Comparison between the compliance and $A^{\mathrm{*}}$ crack measurement methods.

5 CONCLUSIONS

Accurate measurement of crack paths is crucial for validating numerical simulations and improving the reliability of results. In this work, we introduced a novel application of the $A^{*}$ method for tracking cracks in numerical simulations. Unlike conventional methods that only measure the direct path from the crack nucleation site to the crack tip, our approach can track the entire crack path, regardless of its shape. The method is efficient because it searches the nodes closer to the crack path.

We illustrated the effectiveness of the proposed method by presenting fabricated solution cases for both curved and branched crack patterns. Our method was compared to the traditional $A^{*}$ approach, and the results showed excellent qualitative and quantitative agreement. We also presented a time-dependent fabricated case that highlighted the influence of the mesh, crack velocity, and sampling frequency of the results.

Furthermore, we applied our method to crack propagation simulations using the phase-field methodology and obtained an excellent agreement with the compliance technique. This improved accuracy allows for more reliable results and refined parameters.

In future studies, we plan to extend the proposed method to the experimental field using digital image correlation (DIC) methodology, which bridges experimental and numerical analyses.

Acknowledgements:

The authors would like to acknowledge CAPES (Coordenação de Aperfeiçoamento de Pessoal de Nível Superior) for the scholarship under grant 88882.435204/2019-01. J.A. Avila is a Serra Hunter Fellow and a CNPq fellow.

• https://doi.org/10.1590/1679-78257559

References

• Alessi, R., Vidoli, S., and De Lorenzis, L. (2018) A phenomenological approach to fatigue with a variational phase-field model: The one-dimensional case. Engineering Fracture Mechanics, 190:53–73.
• Amor, H., Marigo, J., and Maurini, C. (2009). Regularized formulation of the variational brittle fracture with unilateral contact: Numerical experiments. Journal of the Mechanics and Physics of Solids, 57(8):1209–1229.
• ASTM. E399 17 - Standard Test Method for Linear-Elastic Plane-Strain Fracture Toughness KIc of. Annual Book of ASTM Standards, pages 1–33, 2017.
• ASTM. E647-15e1 - Standard Test Method for Measurement of Fatigue Crack Growth Rates 1. ASTM Book of Standards, 03(July):1–49, 2016.
• Bär, J., (2020) Crack detection and crack length measurement with the DC potential drop method–possibilities, challenges, and new developments. Applied Sciences (Switzerland), 10(23):1–14.
• Bleyer J., Roux-Langlois, C., and Molinari, J.F. (2017) Dynamic crack propagation with a variational phase-field model: limiting speed, crack branching and velocity-toughening mechanisms. International Journal of Fracture, 204(1):79–100.
• Boettinger, W J, Warren, J A, Beckermann, C, and Karma, (2002), A. phase -field simulation of solidification.
• Boldrini, J. L., Barros de Moraes, E. A., Chiarelli, L. R., Fumes, F. G., and Bittencourt, M. L. (2016). A non-isothermal thermodynamically consistent phase field framework for structural damage and fatigue. Computer Methods in Applied Mechanics and Engineering, 312:395–427.
• Carrara, P., Ambati, M., Alessi, R., and De Lorenzis, L. (2020) A framework to model the fatigue behavior of brittle materials based on a variational phase-field approach. Computer Methods in Applied Mechanics and Engineering, 361:112731.
• Creel J.A., Stover S.M., Martin R.B., Fyhrie D.P., Hazelwood S.J., and Gibeling J.C.(2009). Compliance calibration for fracture testing of anisotropic biological materials. Journal of the Mechanical Behavior of Biomedical Materials, 2(5):571–578.
• Dawid, W. and Pokonieczny, K. (2021). Methodology of using terrain passability maps for planning the movement of troops and navigation of unmanned ground vehicles. Sensors, 21:4682, 07.
• Díaz, F. A., Patterson, E. A., Tomlinson R. A., and Yates J. R. (2004) Measuring stress intensity factors during fatigue crack growth using thermoelasticity. Fatigue and Fracture of Engineering Materials and Structures, 27(7):571–583.
• Dijkstra, E. W. (1959). A note on two problems in connexion with graphs. Numerische Mathematik, 1(1):269–271.
• Erke, S., Bin D., Yiming, N., Qi, Z., Liang, X., and Dawei, Z. (2020). An improved A-Star based path planning algorithm for autonomous land vehicles. International Journal of Advanced Robotic Systems, 17(5):1–13.
• Farahani, B.V., Direito F., Sousa P.J., Tavares P.J., Infante V., and Moreira P.P.M.G. (2022), Crack tip monitoring by multiscale optical experimental techniques. International Journal of Fatigue, 155(October 2021):106610.
• Fünfschilling, S., Fett, T., Oberacker, R., Hoffmann, M.J., Ozcoban H., Jelitto, H., Schneider, G.A., and Kruzic, J.J. (2010) R curves from compliance and optical cracklength measurements. Journal of the American Ceramic Society, 93(9):2814–2821.
• Gehri, N., Falcón J.M.-, and Kaufmann W. (2020). Automated crack detection and measurement based on digital image correlation. Construction and Building Materials, 256:119383.
• Hart, P.E., Nilsson, N.J., and Raphael, B. (1968). A formal basis for the heuristic determination of minimum cost paths. IEEE Transactions on Systems Science and Cybernetics, 4(2):100–107.
• Haveroth, G.A., Vale, M.G., Bittencourt, M.L., and Boldrini, J.L.(2020). A non-isothermal thermodynamically consistent phase field model for damage, fracture, and fatigue evolutions in elasto-plastic materials. Computer Methods in Applied Mechanics and Engineering, 364:112962.
• Jobin, T. M., Khaderi S. N., and Ramji M. (2020). Experimental evaluation of the strain intensity factor at the inclusion tip using digital photoelasticity. Optics and Lasers in Engineering, 126(June 2019):105855.
• Kanters, M.J.W., Stolk, J., and Govaert, L.E., (2015). Direct comparison of the compliance method with optical tracking of fatigue crack propagation in polymers. Polymer Testing, 46:98–107.
• Klinsmann, M., Rosato, D., Kamlah, M., and McMeeking, R.M. (2015), An assessment of the phase field formulation for crack growth. Computer Methods in Applied Mechanics and Engineering, 294:313–33.
• Lawande, S.R., Jasmine, G., Anbarasi, J., and Izhar L.I. (2022). A Systematic Review and Analysis of Intelligence-Based Pathfinding Algorithms in the Field of Video Games. Applied Sciences (Switzerland), 12(11).
• Lee, S.Y., Huang, E.W., Woo, W., Yoon, C., Chae, H., and Yoon, S.G., (2015), Dynamic strain evolution around a crack tip under steady- and overloaded-fatigue conditions. Metals, 5(4):2109–2118.
• Liu, J., Liang, X., Xue Y., Fu, Y., Yao, K., and Dou, F. (2020) Investigation on crack initiation and propagation in hydraulic fracturing of bedded shale by hybrid phase-field modeling. Theoretical and Applied Fracture Mechanics, 108(May):102651.
• Martínez-Pañeda, E., Golahmar, A., and Niordson, C.F. (2018). A phase field formulation for hydrogen assisted cracking. Computer Methods in Applied Mechanics and Engineering, 342(March 2020):742–761.
• Mokhtarishirazabad, M., Lopez-Crespo, P., and Zanganeh, M. (2018), Stress intensity factor monitoring under cyclic loading by digital image correlation. Fatigue and Fracture of Engineering Materials and Structures, 41(10):2162–2171.
• Penrose, O. and Fife, P.C. (1990). Thermodynamically consistent models of phasefield type for the kinetics of phase transitions Oliver PENROSE. 43:44–62, 1990.
• Petrini, A.L.E.R., Boldrini J.L., and Bittencourt, M.L. (2021). A thermodynamically consistent phase field framework for anisotropic damage propagation. Latin American Journal of Solids and Structures, 18(1):1–17.
• Saxena, A. and Hudak, S. J., (1978), Review and extension of compliance information for common crack growth specimens. International Journal of Fracture, 14(5):453–468.
• Tortora, G. and Francese, R., (2003). Graph Algorithms.
• Wang, H., Zhou J., Zheng G., and Liang, Y. (2014). Has: Hierarchical a-star algorithm for big map navigation in special areas. In 2014 5th International Conference on Digital Home, pages 222–225.
• Yuan, Y., Ge, Z., Su, X, Guo X., Suo T., Liu Y., and Yu Q. (2021) Crack length measurement using convolutional neural networks and image processing. Sensors, 21(17):1–16.
• Zhang, P., Hu, X., Wang, X., and Yao, W. (2018), An iteration scheme for phase field model for cohesive fracture and its implementation in Abaqus. Engineering Fracture Mechanics, 204(October):268–287.
• Zhang, X., Sloan, S.W., Vignes, C., and Sheng, D. (2017), A modification of the phase-field model for mixed mode crack propagation in rock-like materials. Computer Methods in Applied Mechanics and Engineering, 322:123–136.

Edited by

Editor: Marco L. Bittencourt and Josué Labaki

Publication Dates

• Publication in this collection
  17 July 2023
• Date of issue
  2023

History

• Received
  15 Mar 2023
• Reviewed
  24 Mar 2023
• Accepted
  30 May 2023
• Published
  05 June 2023