ANSYS implementation of the phase field fracture approach

In this study, we present a new implementation of the phase field fracture approach in the finite element code ANSYS and its numerical background. The framework is general, and is supported by addressing several classical 2D boundary value problems as well as the ductile fracture and 3D surface flaws behaviors of particular interest. The 3D implementation exploits the analogy between the phase field formulations and the magnetic vector potential equation. The influence of the mode mixity and biaxiality loading conditions of the cracked bodies on phase fields is evaluated as a function of the crack length scale parameter, characterising the scale at which damage effects become significant. As a result of the FE calculations of phase field distributions for propagating cracks, the effects of both the fracture mode and the biaxial stress-strain state are determined. The size of the fracture process zone or damaged region is determined across a wide range of cracked body loading conditions. Developed code: https://github.com/Andrey-Fog/ANSYS-USERELEMENT-PHFLD.


INTRODUCTION
he phase field fracture approach in the recently literature enjoys great popularity and has been successfully applied to model brittle fracture, ductile damage, hydrogen assisted cracking, fatigue crack propagation, to name a few.This concept has enabled predicting complex fracture phenomena such as crack initiation, coalescence, branching, kinking and bifurcation.Phase field methods for brittle fracture builds upon from the principal study of Francfort and Marigo [1], who were the first to employ the Griffith's thermodynamic framework into variational formulations by treated elastic fracture as an energy minimization problem.According to the Griffith's energy balance crack growth will take place if a critical energy release rate is attained.Bourdin et al. [2] regularized the discrete crack topology by means of a scalar damage variable and a diffuse crack representation.This variable is termed as the phase field, or phase field order parameter.In this, brittle fracture is numerically treated as a coupled, i.e., displacement and phase field problem.By Miehe et al. [3] has been introduced a phase field law primarily to model brittle fracture which than was employed in several studies to model ductile fracture.

T
The variational formulation was further modified and extended to multi-dimensional mixed-mode brittle fractures also targeting the response of visco-elastic solids and ductile materials.The phase field representation of fracture successfully has been extended to the ductile regime also within the context finite strains and the employment of different constitutive equations describing the nonlinear material behavior.By the authors Ambati et al. [4], Borden et al. [5] several methods and algorithms of numerical implementation have been developed to tackle to model fracture of ductile materials.Very recently, the finite element-based phase field damage method formulations have been introduced within a coupled strain gradient plasticity problem and hydrogen assisted cracking framework.Based on a new non-conventional thermodynamics approach, phase field method for modelling ductile fracture in damaged solids has also been established by Tsakmakis and Vormwald [6].In addition, phase field models have been applied successfully in fracture of visco-elastic and ductile materials [7,8].The application of phase field fracture approaches has also triggered a notable interest to solve the coupled deformationdiffusion-fracture problem [9].One of the obvious and expected applications of the theory of phase field fracture was the description of the fatigue crack propagation.In this case the damage resulting from the application of cyclic loads is modelled by means of a fatigue degradation function which initially has been introduced in work [10].Subsequently, applications to hydrogen-assisted fatigue [11], shape memory alloys [12], crack growth [13] and general fatigue problems [14] were demonstrated.Along with general formulations and numerical justification the phase field method has found application in the simulation of fractures in 3D problems [5].In Wu et al. study [15], several 3D benchmark problems involving mixed mode I/II/III failure in brittle and quasi-brittle solids has been addressed.Yin and Kaliske, Yin et al. have shown the performance of the proposed phase-field formulations of fracture by means of representative 3D numerical examples [7].Interest in the phase field fracture method has increased significantly in recent years, and comprehensive reviews on this issue have already been published by Wu et al. [16].Shlyannikov et al. [17] supplemented this review with a discussion about the variability and behavior of the experimentally determined critical energy release rate or fracture toughness of the material G c under monotonic and cyclic loading.Attempts to experimentally validate the phase field fracture method have also been provided [6].By the authors [18] in order to study fatigue mechanisms numerically, a simple strategy to incorporate residual stresses in the phase-field fatigue model is presented and tested with experiments.To introduce the fatigue effects, Carrara et al. [10] propose to modify the fracture energy by introducing t as the pseudo-time, α(t) as a properly defined cumulated history variable, whereas the history variable f(α(t)) is a fatigue degradation function.In this framework, the dissipated energy is a process-dependent quantity and is no longer a state function.This allows the inclusion of fatigue effects as a reduction in the local fracture toughness of the material proportional to the cumulative history variable α(t).The fatigue degradation function f(α(t)) describes how fatigue effectively reduces the fracture resistance of the material G C [10].The fatigue-related quantity α(t) is a threshold parameter controlling when the fatigue effect is triggered and it is meant to be material parameter to be determined on an experimental basis.To this end it is possible to apply the peak strain εy as functions of the Young's modulus E, the regularisation length ℓ and the fracture toughness G C [19,16].Using these parameters leads to a direct dependence of the threshold value of the fatigue degradation function on the fracture toughness G C .Initially, the well-known commercial finite element method code ABAQUS was the main platform for implementing phase field fracture technologies [20].The phase field is solved for at the finite element nodes, as an additional degree of freedom.Martínez-Pañeda and Golahmar provided an efficient and robust implementation of the phase field method in the commercial finite element package Abaqus, enabling to model interactions and branching of cracks of arbitrary topological complexity [21].Later Navidtehrani et al. has been presented a new ABAQUS implementation [22].The aim of this work was to present, establish the features and background the implementation on a wide range of diverse numerical examples of the phase field method in the ANSIS computer code environment.In this paper a new description is proposed for an implementation of the phase field fracture formulations into the ANSYS finite element software, which differs from the commercial FEM package ABAQUS.This manuscript is structured as follows.Section "Theoretical background" provides a concise overview of phase field fracture theory.Subsequent sections include "Finite element implementation" in ANSYS and "Representative results".Initially, we address fracture using a benchmark example of a cracked body under uniaxial tension.Then, we analyze crack growth conditions with a cracked plate subjected to shear.Thirdly, the compact tension (CT) specimen behavior is analyzed, and is shown the results of 2D solutions by various elastic and ductile constitutive equations choices.Than the compact tension shear (CTS) specimen is considered in order to assessment of mode mixity effect.Finally, we simulated the 3D through-thickness crack and surface flaw propagation, which is very important for structural integrity assessment of the advanced materials.A comprehensive 3D analysis is performed, including the modelling of an inclined surface crack subjected to biaxial loading as well as interaction and coalescence of two semi-elliptical flaws.The manuscript ends with concluding remarks in Section "Conclusions".

THEORETICAL BACKGROUND
he computational simulation of solid material failure due to fractures characterized by sharp crack discontinuities encounters challenges when dealing with complex crack geometries.These challenges are effectively addressed through the use of a diffusive crack modeling approach that incorporates a crack phase-field.Central to these models is the approximation of sharp crack profiles by diffuse cracks spread over a finite-width localization band controlled by the length scale l (Fig. 1).Researchers such as Francfort, Marigo, and Bourdin et al. have introduced dissipation functions linked to the evolution of the crack surface functional, which is assumed to exhibit complete dissipative behavior [1,2].These dissipation functions are formulated to depend on the phase-field rate and its gradient, embodying a gradient-type dissipative material model.They are designed specifically to simulate crack propagation by facilitating localized phase-field growth.The phase field paradigm is an auxiliary (phase) field variable  that was defined to describe discrete discontinuous phenomena, such as cracks, with a smooth function.The key aspect of PFF models is the approximation of the sharp crack topology by a diffuse crack smeared within a localisation band of a finite width controlled by the length scale l.The key idea is to smear a sharp interface into a diffuse region using this phase field order parameter , which takes a distinct value for each of the two phases (e.g., 0 and 1) and exhibits a smooth change between these values near the interface.From the viewpoint of material modeling, a phase-field approach to fracture is conceptually in line with models of continuum damage mechanics (CDM), a discipline pioneered by Kachanov, where a scalar damage field may be interpreted as the phase-field.In this sense, the subsequent framework of regularized brittle crack propagation may be considered as a gradient-type damage model with particular definitions of an energy function with a gradient-type regularized surface energy.Accordingly, the loss of stiffness associated with mechanical degradation of the material is characterized as a function of .This is done by means of the so-called degradation function g(), which relates the stored bulk energy per unit volume to the strain energy density of the undamaged solid.Clearly, such an approach is related to the continuum damage mechanics, where the regularising length scale parameter l dictating the size of the fracture process zone or the size of the damaged region.This regularisation parameter is a material property governing the material strength and characterising the appearance and growth of micro-cracks and micro-voids.With this viewpoint in mind, we performed a special study using a Merlin Zeiss backscattered scanning electron microscope (SEM) to identify the fracture process zone and the phase field on the outer surface of the tested high-strength steel specimen.

T
In Fig. 2 backscattered SEM images of fracture surfaces of specimen are displayed which is clear evidence of phase field paradigm.The main crack in Fig. 2a is surrounded by several smaller additional cracks, the distribution density of which decreases with distance from the main crack.In addition, Fig. 2a is a good image of the fracture process zone in the vicinity of the main crack.It should also be noted that at high magnification the crack tip has a finite radius of curvature (Fig. 2b).These findings support the damage model based on the concept of the phase field theory concerning a diffuse crack smeared within some finite area controlled by the length scale factor l.
Within the phase field modeling framework, the interface is represented over a diffuse region using an auxiliary field variable  [3].This variable adopts distinct values for each phase  = 0 indicating no crack and  = 1 representing a fully fractured state at the integration point), transitioning smoothly between these values near the interface.In definition (x) is the crack phase field with Dirac delta function.The nonlocal damage variable  near the crack surface is approximated by the exponential function (1): ( ) e x l x   Eqn.
(1) depicts a regularized or diffuse crack topology (Fig. 1), where l is the length scale parameter controlling the width of the diffuse crack zone.It is crucial to recognize that l does not directly equate to the actual crack width, as the crack is distributed smoothly across the entire domain.The subsequent discussion will reveal that the phase field variable  evolves over time according to a partial differential equation (PDE).This methodology allows for the simulation of complex interface evolution phenomena by solving a system-wide set of PDEs, thereby removing the necessity for explicit interface condition management.

Representation of crack surface density function
The crack surface is related to the crack length scale parameter l (Fig. 3).The crack surface density is introduced by means of the regularized crack functional as [21]: where (, ) is the crack surface density function in perpendicular to the crack surface direction.In higher dimensions it is expressed as: Figure 3: The diffusive crack surface  l () is a functional of the crack phase field.

Governing balance equations of coupled problem
The total potential energy functional of a solid body is presented in the following form: where  b (u, ) -is the stored bulk energy, and  s () -is the fracture energy dissipated by the formation of a crack.The total potential energy functional adopts the form incorporating the parabolic degradation function g()=(1-) 2 , which characterizes the decay of stored bulk energy due to damage evolution: where () is the strain energy density.

The total potential energy functional for ductile fracture problem
The total potential energy functional of a solid body is presented in the following form: and it depends on the displacement u and the fracture phase field .The energy storage function describes the stored bulk energy of the solid per unit volume.
where  0 is the strain energy density, g()=(1-) 2 is the degradation function.The total strain energy density includes both elastic  el and plastic  pl parts: According to the principle of superposition the elastic part in (8) can be easily obtained by summing for all directions: To impede the formation and propagation of cracks under compressive stress, the strain energy density can be partitioned into separate tensile and compressive components in the following manner: Preventing damage under compression involves decomposing the strain energy density, typically addressed through two distinct methods.One such method is the volumetric-deviatoric split, introduced by Amor et al. [23], formulated as follows: where  el represents the deviatoric part of the elastic strain tensor,  and  Lame parameters, tr -denotes the trace of the matrix, and ⟨ ⟩ signifies the Macaulay brackets, K -bulk modulus.The second one is the so-called spectral decomposition, proposed by Miehe [3], The second approach is the spectral decomposition, introduced by Miehe [3], which utilizes the spectral decomposition of the strain tensor   =, ⟨ I ⟩  n I  n I with  I and n I being, respectively, the principal strains and principal strain directions (with I = 1, 2, 3).The decomposition of strain energy is then expressed as: Damage is an irreversible process:   0. To enforce irreversibility, a history field variable H is introduced, which must satisfy the Karush-Kuhn-Tucker (KKT) conditions: 0, 0, ( ) 0 Accordingly, for a current time t, over a total time , the history field can be defined as, To determine the plastic part, it is necessary to introduce a yield function.The Voce exponential law was taken as the yield function: where To calculate the strain energy density, it is necessary to record the obtained values of the strain energy density, strain, and stress at the previous time increment t.The strain energy density is updated with the current time step t + ∆t, Thus, the plastic part of the strain energy density: Consider the minimization problem of the internal potential energy with respect to the phase field.It yields the strong form of the evolution of the crack phase field along with the Neumann-type boundary condition:

Numerical implementation 2D problem in ANSYS
o implement the phase field model in the finite element method, it is necessary to introduce an additional degree of freedom.The phase field model is implemented by means of Ansys UEL subroutine which allows for user-defined computation of the element tangent stiffness matrices and the nodal force vectors.We consider isoparametric 2D quadrilateral elements (linear and quadratic) with 3 degrees of freedom per node, i.e. u, v and , and four integration points.
To solve the problem using the finite element method, it is necessary to obtain a system of linear ones in the following form: where K -stiffness matrix of the element, {u} -degrees of freedom vector (displacement field by default) and {r} -loads vector.The simplest way to implement a phase field model is to replace the degree of freedom u with  for 2D case.
T Using Voigt-notation in a 2D space, the displacement field u={u, v} T and the phase field  in ( 5) can be discretized as In ( 20) N -denotes the shape function associated with node i from all m nodes of element.The corresponding derivatives can be discretized as: where  = { xx ,  yy ,  xy } T .The strain-displacement matrices are expressed as: where N i,x and N i,y are the derivatives of the corresponding shape function with respect to x and y, respectively.According to ANSYS documentation [24], the displacement vector in ( 19) must be present in the following form {u}={u 1 ,v 1 , 1, … , u 4 ,v 4 , 4 } , and the nodal load vector {r}={r u1 ,r v1 ,r 1, … , r u4 ,r v4 ,r 4 } .Thus, for implementation it is necessary to arrange the components of the stiffness matrix according to: The elements of stiffness matrix (23) and load vector can be calculated from: for degrees of freedom u i and v i , and for degrees of freedom  i .
Employing the Newton-Raphson method, the final system of equations can be solved iteratively: The use of a Newton-Raphson backward Euler scheme implies that the numerical solution is unconditionally stable, i.e., achieving convergence implies attaining the equilibrium solution.The extension to a three dimensional case is straightforward.

ANSYS particularities of 3D problem for degrees of freedom
For 3D problems FEM code ANSYS does not allow a free choice of degrees of freedom.Therefore, only built-in sets should be used.According to Navidtehrani et al. [22], the implementation utilized an analogy between the evolution law of the phase field and the magnetic vector potential.This approach enables using the vast majority of ANSYS' in-built features, which help coding the user-defined coupled field elements.
The following degrees of freedom can be used for this purpose in ANSYS finite element library for DOF command: UX, UY, UZ (structural displacements); TEMP (temperatures); PRES (pressure); VOLT (voltage); MAG (magnetic scalar potential); AZ (magnetic vector potential); CURR (current); EMF (electromotive force drop); CONC (concentration) and more.Any of these DOD can be used as a phase field degree of freedom.The basic equations of the vector potential formulation of the finite element method are: where {v m } -magnetic vector potentials at the nodes of the element (input/output as AZ), J i -current density vector (input/output as CSGZ), [ C ] is damping matrix and [ K ] is coefficient matrices.The analogy of this Eqn.(27) with the phase field evolution law ( 19) is evident, with the magnetic vector potential acting as the phase field {v m }=.Assuming that the final coupled system of equations at each integration point for 3D case can be reformulate for the phase field local force balance in the following form: where final stiffness matrix and load vector can be calculated from ( 24) and (25).Eqn. ( 29) is the basis of the proposed user element for 3D phase field fracture problems which can be solved by a default ANSYS solver.In the present study this elaborated method was implemented in ANSYS software via a user programmable feature to solve the represented below 3D problems.As was mentioned above ANSYS allows selecting different degrees of freedom for user programmed element.It should be noted that the default convergence criteria do not lead to a satisfactory solution.To set the correct convergence parameters for nonlinear analyses the CNVTOL command should be used.The Tab. 1 correlates the list of the degrees of freedom with the nodal load vectors that can be used for phase field fracture problems without changing of source code of the userdefined element.

REPRESENTATIVE RESULTS
he effectiveness and capabilities of the current ANSYS implementation will be exemplified through simulations of fracture in several foundational boundary value problems.Initially, we will explore the initiation and propagation of cracks in a notched square plate under both uniaxial tension (as outlined in the subsection "Notched square plate subjected to tension test") and shear loading conditions (as detailed in the subsection "Single Edge-Notched Shear Test").Then, the failure of the compact specimen (CT) subjected to tension, under elastic and elastic-plastic deformation, is simulated in subsection "Compact tension specimen".Next subsection "Compact tension shear specimen" is devoted to analysis of the behavior of CTS sample under mixed mode loading.Finally, in subsection "Full 3D problem", 3D models of a notched square plate and an inclined surface crack subjected to biaxial loading are developed, including the two semi-elliptical flaws, to determine the nucleation and coalescence of cracks.Each of the listed examples will be accompanied by a table containing a list of parameters used in the corresponding numerical calculations.This list includes the following parameters that are contained in the main constitutive equations (1,3,4,13,15) used: E is Young's modulus;  is Poisson's ratio; l is the phase field length scale; G c is the material fracture toughness;  0 is the yield stress; R inf is exponential coefficient; b is exponential saturation parameter;  is crack orientation angle; H min is minimum value of the fracture driving force; h is the characteristic element size; u is a displacement-driven deformation by prescribed incremental displacements; a/c is aspect ratio of semi-elliptical crack.

Notched square plate subjected to tension test (Single Edge-Notched Tension Test)
Initially, we will explore the scenario of unstable crack growth in a square plate with a notch subjected to uniaxial tension.This case is widely recognized as a classic benchmark in the study of phase field fracture mechanics, as established in the pioneering works by Miehe et al. and Martínez-Pañeda and Golahmar [3,21].The phase field model is implemented via the ANSYS USERELEM subroutine, which enables user-defined calculations for the element tangent stiffness matrices and the nodal force vectors.The geometric configuration and boundary conditions are depicted in Fig. 4a.The specimen is specifically subjected to mode I fracture conditions, characterized by a prescribed vertical displacement at the remote boundary.The specimen is subjected to a displacement-controlled deformation through incremental displacements as it approaches the peak load.The parameters used in the analysis are detailed in Tab. 2. We discretise the model using linear isoparametric 2D quadrilateral elements with 3 degrees of freedom per node, i.e. u, v and , and four integration points.In Fig. 4b, the mesh is tailored to follow the anticipated crack path, ensuring that the characteristic element size is at least four times smaller than the phase field length scale l.This approach was taken to T guarantee mesh independence in the region where crack propagation occurs.Approximately 10,000 quadrilateral elements were used in total.
(a) (b) The force-displacement relationship for a notched square plate under tension testing is shown in Fig. 6a, contrasting elastic and ductile material behaviors.This comparison highlights how the phase-field formulation's choice of constitutive equations significantly influences peak load outcomes.The fracture patterns at different loading stages are visualized in Fig. 6(b,c).The response is not symmetric, as the lower bound is fully clamped, and the crack is rather diffuse, as expected given the choice of governing equations.In addition, the phase field damage is extending at an angle of approximately 40-45 degrees relative to the plane of symmetry of the sample.The force-displacement relationship for a notched square plate under tension testing is shown in Fig. 6a, contrasting elastic and ductile material behaviors.This comparison highlights how the phase-field formulation's choice of constitutive equations significantly influences peak load outcomes.The fracture patterns at different loading stages are visualized in Fig. 6(b,c).In order to point out the effects that arise due to the length-scale parameter l and the constant ratio l /h, different simulations with elastic model are performed.For fixed critical energy release rate G c = 2.7 MPa•mm, the influence of the length-scale parameters is analyzed.The predicted crack path is showcased in Fig. 8 by plotting the contours of the phase field variable .The phase field response of the crack is rather diffuse, as the absolute value of the length-scale parameter is increased, as expected given the choice of l from 0.02 to 0.12.It can be seen that, the sharpest crack pattern is obtained for the smallest length-scale parameter l = 0.02mm.The subsequent study analyzes the influence of the discretization on the global response.
To verify objectivity, the simulation discussed in Fig. 5 is repeated with a finer discretization of 20,000 elements, resulting in an effective element size of approximately h ≈9×10 -3 mm.Comparing both structural responses confirms that the results are independent of the mesh size.

Single Edge-Notched Shear Test (Crack branching under mixed mode fracture)
Crack branching remains one of the most complex issues in brittle fracture mechanics.This phenomenon has been previously examined by Bourdin et al. [2] and Tsakmakis and Vormwald [6].The literature suggests that a bifurcation occurs at the crack tip, leading to competing paths for crack propagation.Some researchers assert that the crack tends to deflect towards the bottom-right corner of the sample, while others propose that two branches form, inclined relative to the horizontal axis.The potential directions of crack propagation have been linked to several factors, such as the anisotropic damage evolution within the crack tip fracture process zone, the plastic characteristics of the material, and the specific phase field models employed.These aspects remain open questions in phase field fracture for both elastic and ductile materials, as they lack extensive experimental validation.
In the subsequent investigation, we aim to analyze stable crack growth through the simulation of fracture in a notched square plate under shear loading.This study builds upon the framework discussed in the subsection "Notched square plate subjected to tension test" The geometric configuration and boundary conditions are depicted in Fig. 9. Specifically, a horizontal displacement is applied at the top edge of the plate, while the bottom edge is fully constrained.The crack propagates parallel to both the upper and lower edges of the plate, with fixed orientations prescribed along these boundaries.The initial dimensions of the crack and the size of the specimen remain unchanged, consistent with those utilized in the previous study on uniaxial tension.A discretization of 20,000 isoparametric 2D quadrilateral elements is used with mesh refinement along the expected crack path.The initial data used in this calculation are given in Tab. 4. It should be noted that the material fracture toughness G C depends on the stress state, which can be expressed in different values of G CI , G CII , and G CIII , or as a generalized dependence on the phase angle, G C(α).In the case of considering specific problems of pure mode I or mixed mode fracture, the values of the equivalent parameter as G eqv = f(G I , G II , G III ), which is a function of all three fracture modes, should be used to obtain comparative quantitative estimates of the behavior of the material for each specific situation.The dissipation function is governed by the magnitude of the geometric crack function and, consequently, a scaling constant l = 0.04 mm.For the AT1 model, which was used in the calculations () = , c  =2/3 and minimum value of the fracture driving force   The load-displacement curves for two branches and one inclined cracks are depicted in Fig. 10.Despite the differences in crack growth trajectories, the behavior of the load-displacement curves is close to each other and shows moderate descent of the load in the post-critical regime with respect to mode I elastic solution (Fig. 8a).The fracture patterns at several loading stages are visualized in Fig. 11 for single-edge-notched plate subjected to shear test and are in good agreement with the results reported in the literature [2,4,6].A comparison of phase field fracture fields as a function of applied force direction indicates that the deviation of the crack angle from  =0˚ to  =7˚ has sufficient effect on the damage fields or crack trajectory.It should be noted that the state of initial pure shear in fracture mechanics is characterized by extreme instability in the sense of the direction of possible crack propagation.In the case of  =0˚, the phase field response in Fig. 11(a-b) exhibit two symmetric branches of a preexisting crack, whereas at  =7˚, crack kinking is observed, typical for a mixed mode fracture as can be seen in Fig. 11(c-d).
In the next example, the influence of the elastic-plastic material behavior on the predicted crack path is examined with respect to linear elastic stress strain state.The material response in this simulation is supposed to exhibit the isotropic hardening according to Eq.( 15) and Tab. 4. The examined example of a single-edge-notched plate subjected to two types of loading angles  =90˚ and  =0˚, illustrates the material behavior under classical modes of failure in fracture mechanics, specifically pure mode I and pure mode II.Fig. 12 shows a comparison of the load-displacement curves for these loading conditions as a result of solving elastic (Fig. 12a) and plastic (Fig. 12b) problems.From these data it follows that shear mode II leads to a decrease in the peak load and an increase in accumulated strains compared to mode I.
(a) (b) In the literature on phase fields fracture, there are a number of classical benchmark examples that are used to substantiate the formulation of new theories and basic equations, the robustness and capabilities of the proposed algorithms, as well as numerical methods in their implementation.Traditionally, simplified models, configurations, and boundary conditions are considered for these purposes.Clearly, caution must be exercised in drawing conclusions from these modelling solutions in terms of practical applications.The results obtained from parametric studies are important as useful recommendations rather than results of direct practical application.In this regard, the two examples discussed below in this work relate to specific sample geometries and their loading conditions, which are widely used in experimental fracture mechanics.

Compact tension specimen (plane problem)
The first example which is addressed to practical applications is a compact tension (CT) sample.The geometry is shown in Fig. 13, with the dimensions given in mm.The load is applied by prescribing the vertical displacement of the nodes in the pin holes, which are not allowed to move in the horizontal direction.We assume a material with main mechanical properties listed in Tab. 5.The objective now is to evaluate the impact of different constitutive models on the phase field fracture description.As outlined in the subsection "The total potential energy functional for ductile fracture problem," we consider two options: elastic and ductile models.We begin by comparing the predicted fracture patterns in a CT specimen and a notched square plate under uniaxial tension (subsection "Notched square plate subjected to tension test").In Fig. 14a, force-displacement curves are plotted for the CT specimen using two sets of material properties: elastic and elastic-plastic states.As observed, damage leads to a moderate decrease in load with significant deformation across the specimen section.However, predictions for the singleedge-notched plate under tension exhibit a different trend (Fig. 6a).Specifically, the load at crack initiation appears highly sensitive to the choice of constitutive equations, with lower strength values associated with the ductile model.This sensitivity arises due to the influence of boundary conditions on phase field fracture, particularly when the lower edge of a notched plate is fully clamped.The resulting fracture patterns for compact tension specimen at different loading steps are illustrated in Fig. 14(b,c).It should be noted that the phase fields of fracture in the CT sample under consideration are the same in type for elastic and plastic deformation.In both scenarios, the qualitative results align well with extensive experimental data.As depicted in Fig. 14(b,c), the response exhibits symmetry, and the crack propagates steadily along the anticipated mode I path.In contrast, when using the same governing equations, the model problem of single-edge-notched plate under tension during plastic deformation shows a deviation of the crack path from the plane of symmetry (Fig. 6c).Obviously, this is due to the different configuration and method of applying the load for the plate and the experimental sample.

Compact tension shear specimen
The impact of different phase field models and loading conditions discussed earlier is explored through modeling crack propagation in a Compact Tension Shear (CTS) sample.The dimensions of the sample, given in millimeters, are depicted in Fig. 15 along with the utilized test setup.S-shaped grips are utilised to achieve mode mixity between pure mode I and II loading conditions (Fig. 16a) together with the finite element mesh (Fig. 16b).Pure mode I is obtained when force F 1 is applied in the  = 90˚ direction, whereas pure mode II is achieved by applying force F 2 in the  = 0˚ direction.The specimen is subjected to symmetric displacement at the loading pins.The contact between the pins, S-shaped grips and the specimen in all eight holes is defined as surface to surface contact with a finite sliding formulation.For tensile and shear test conditions loads are applied by prescribing the vertical u y = 0.15 mm and horizontal u x = 1 mm displacements through the pin holes mixed mode grips, respectively.A total of 19,521 4-node plane strain elements are used, with the characteristic element size in the crack propagation region being 5 times smaller than the phase field length scale.The assumed values of a material elastic properties and loading conditions are given in Tab. 6.We consider isoparametric 2D quadrilateral elements quadratic with 3 degrees of freedom per node, i.e. u,v and , and four integration points.A total of 50000 quadrangular elements are used for the tensile test and 200000 for the pure shear loading conditions.The situation of initial pure shear is a subject of ongoing discussion in the research community.In continuation of this discussion, this paper presents an experimental substantiation of our numerical results.In the literature, the main question is which damage mechanism is dominant -branching or kinked.Following the formulation for the total potential energy functional (subsection "The total potential energy functional for ductile fracture problem"), damage can be reads by decomposition of the strain energy density as volumetric-deviatoric split (Eq.11) or spectral decomposition (Eq.12).Fig. 17 shows the experimental trajectories of crack growth obtained by Shlyannikov and Fedotova in the CTS sample for the conditions of pure mode I (Fig. 17a) and the initial mode II of pure shear (Fig. 17c, f) [25].These results relate to samples made of high-strength steel 34X.It can be noted that for the same loading conditions and elastic properties of the material, branching and kinking of the initial pure shear crack are possible (Fig. 17c,e).The results of modeling phase fields for the corresponding loading conditions are illustrated in Fig. 17(b,d,f).The phase fields in Fig. 17e correspond to the use of Eq.9, whereas in Fig. 17f the phase fields refer to Eq.12.Recall that the differences in these equations lie in the decomposition of the strain energy density.Thus, due to various forms of writing the components of the governing equations of the energy balance, it is possible to obtain an adequate description of the experimental data after the fact, however, the question of predicting the crack growth trajectory remains open.In the context of the following elastic calculations, it is essential to establish the most physically meaningful relationship between G c , and l, considering that the length scale dictates the size of the fracture process zone and therefore cannot be regarded merely as a numerical parameter.A similar conclusion about the significance of the choice of constitutive equations was also expressed by Tsakmakis and Vormwald [6] based on a comparison of isotropic and kinematic hardening models in the analysis of phase fields during nonlinear cyclic deformation.
(a) ( The loading-displacement evolution regarding the phase-field models according to Eq.9 and Eq.12 are shown in Fig. 18 at three different fracture status.It can be clearly observed that the loading-displacement relation is a function of mode mixity showing higher peak loads for initial pure shear mode II ( = 0˚) by employing of Eq.( 9) compared to the pure mode I tension condition ( = 90˚).In addition, the slope angle of the curves in the first deformation section before reaching the peak load is also different, with a smaller value during shear (Fig. 18a) for branching cracks.It is obvious that these states correspond to elastic modulus under tension and pure shear.Noteworthy is the different behavior of the curves during branching and kinking of cracks, the trajectories of which relate to mixed mode fracture.It can be seen in Fig. 18b that the profiles of the load-displacement curves regarding such mode mixity also show differences.An abrupt failure is observed for both cases at the loading direction  = 0˚, nevertheless, the specimen with branching dominant damage mechanism according to Eq.( 9) yields a smaller value of peak load with respect to kinking crack path behavior in CTS specimen describing by Eq. (12).Comparing these simulation results to Shlyannikov and Fedotova experimental data both the predicted crack paths and the profiles of the loading-displacement relations show thorough agreement [25].

Full 3D problem
Now we showcase the potential of the framework presented demonstrating the behavior of phase field fracture in 3D solids.We do so by simulating the fracture of brittle solids with through-thickness and surface crack subjected to uniaxial tension or biaxial tension-compression loading conditions.The results presented below were obtained in the order of implementing the developed algorithm in the ANSYS software, taking into account the features of solving three-dimensional problems, which are described in subsection "ANSYS particularities of 3D problem for degrees of freedom".Next, we investigate the mode-I fracture problem through a three-dimensional tension test.The geometric configuration and boundary conditions are shown in Fig. 19(a).The bottom edge is fixed, while vertical displacement is applied to the top edge.A notched plane specimen is discretized using an unstructured mesh of approximately 100,000 quadrilateral elements.
We utilize isoparametric 3D-Brick elements that are quadratic with 4 degrees of freedom per node (u, v, w, and ) and feature four integration points.The mesh is particularly refined in the region where crack propagation is anticipated (Fig. 19(b)).
(a) (b) (c) As expected, the crack starts at the notch tip and propagates straight through the specimen.The crack shows some roughness but is almost planar.
In the next example we simulate a semi-elliptical inclined surface crack in a plate subjected to a remote uniform biaxial loading (Fig. 21a) for a linear elastic material.This 3D example demonstrates that phase-field models have the potential to predict complex failure behavior of three-dimensional structures under severe load conditions.The biaxial load ratio,  =  xx / yy , defines the relationship between the nominal stresses.In Fig. 21b, the plate thickness featuring the inclined surface flaw is uniformly set to t=0.25 mm across all examined scenarios.The surface flaw depicted in Fig. 21c experiences mixedmode fracture conditions, characterized by the presence of all three stress intensity factor components: mode I, mode II, and mode III.This subsection primarily aims to examine the effect of the aspect ratio  elp = a/c of a semi-elliptical surface flaw and the biaxial loading ratio  on phase field behavior, with an initial crack angle consistently set at  = 45˚.
(a) (b) (c) Figure 21: A semi-elliptical inclined surface crack in a plate under remote uniform biaxial loading.
In the phase-field method, crack evolution is a result on the global level, that is based on the formulation and evaluation of the local strain energy density, driving force and the material fracture resistance.Therefore, the exact quantity of strain energy density driving the crack propagation must be influenced by both the strain state and the orientation of the crack surface in relation to the remote uniform biaxial load.The combination of the spatial description of the crack surface with the governing physical principle of crack surface formation, i.e. strain energy dissipation, provides an approach, that can model and predict the propagation of cracks along realistic and experimentally validated paths.In the present study, several 3D benchmark problems involving mode I, I+II or I+III failure in brittle and quasi-brittle solids are addressed based on our recent numerical [26] and experimental [27] progresses.
(a) (b) We proceed now to simulate the fracture of a plate subjected to biaxial tension or tension-compression.Fig. 22 shows the geometry of a plate with an inclined surface crack and boundary conditions for implementing biaxial loading.The crack was modeled in such a way that the aspect ratio of the initial crack was a/c = 1 and a/c = 0.5, where a = 0.025 mm, c = {0.025,0.05}mm.Three different cases for combining biaxial loading and the initial orientation of the surface of a semi-elliptical crack are considered.The first case represents pure mode I under equibiaxial tension when  = +1 and  = 45˚.The second situation refers to pure mode II under equibiaxial tension-compression with  = -1 and  = 45˚.The third case is a typical mixed mode fracture under uniaxial tension  = 0 of the initial surface crack located at an angle  = 45˚.The material properties are taken to the computations listed in Tab. 7.  It should be noted that in the three-dimensional setting, the number of elements grows dramatically therefore a simplified model has been used for the analysis as shown in Fig. 22b.The numerical model is uniformly discretized by 160000 quadrilateral 8-node linear brick elements.The 3D plate is subject to a displacement-driven deformation by prescribed incremental displacements u.
(a) (b)  [27].As anticipated, a plate initially containing a surface defect exhibits greater loadbearing capacity under uniaxial tension compared to biaxial loading conditions.Specifically, for plates with pre-existing surface defects, the stiffness of the material degrades more rapidly when the defect is semi-elliptical rather than semi-circular, leading to differences in the maximum force attained in each case.In addition, mixed modes of failure at  = 0 show a higher peak load compared to the pure tensile mode I at  = +1 and the initial pure shear mode II at  = -1.The resulting force versus displacement response reveals a quantitative agreement with the computations by Shlyannikov and Tumanov [26].
The surface crack growth trajectories predicted for the three biaxial loading cases described above are shown in Figs.Under uniaxial tension at  = 0 of a plate with an inclined crack  = 45˚ (Fig. 25), the dominant mechanism is mixed modes of fracture, in which the direction of crack growth changes along a curved semi-elliptical front.It can be noted that, first of all, the phase fields of damage propagate in an inclined plane reaching the edge of the plate and after that they cross the entire thickness and the defect becomes through.In this case, the contours of the phase fields have direct symmetry in both the horizontal and vertical planes and show reverse symmetry on the free surface of the plate.
Under conditions of equi-biaxial tension-compression at  = -1 and  = 45˚ (Fig. 26), the initial semi-elliptical defect demonstrates the dominant mechanism of out-of-plane shear.The phase fields in Fig. 27 show that the surface crack during its growth at the initial stage quickly becomes through-thickness with small sizes on the surface of the plate.Further propagation of the crack occurs as a through-thickness defect, predominantly in the horizontal plane.In this case, the fracture patterns of the phase fields in all mutually perpendicular planes are symmetrical relative to the coordinate axes.The obtained numerical phase fields confirm recent experimental data [27], which indicate that surface defects grow at different rates on the free surface of the sample and at the deepest point of the semi-elliptical crack front.The phase field contours shown in Figs.24-26 present a singular perspective for scrutinizing the progressive accumulation of damage along the fronts of semi-elliptical cracks, contingent upon diverse biaxial loading conditions.Three-dimensional effects most clearly affect the behavior of the mode mixity parameters, which convincingly prove, within the same crack front at  = -1 and  = 45˚, a change in fracture modes from pure mode II at  = 0˚on the surface of the sample to pure mode III at  = 90˚the deepest point of the front.In this way, the processes of defect evolution from a part-through-thickness surface crack to a through-thickness crack, which completely intersects the section toward the thickness direction of the sample, are visualized.
Structural integrity assessment often involves estimating surface flaw growth and coalescence from initial operation damages as a function of time.One major merit of phase-field models for fracture is that cracks nucleation, propagation, branching, merging, coalescence and even fragmentation can be accounted for within a standalone regularized variational framework.We investigate the capabilities of the present phase field formulation to predict complex crack propagation due to existing defects.To demonstrate the capability of the phase field approach implementation in ANSYS code, we now consider a 3D example in which two cracks interact and coalesce.We consider a plate with two preexisting semi-circular cracks, as outlined in Fig. 27a, which placed in horizontal section at middle height and symmetrically with respect to the vertical axis of the specimen.The geometry and boundary conditions are illustrated in Fig. 27a.The thickness of the plate is taken as t =0.3 mm.Two initial cracks were modeled in such a way that the aspect ratio of the initial crack was a/c = 1 and the semi-axes are a = c = 0.05 mm.We assume E = 210 GPa,  = 0.3 and energy release rate G c = 2.7 MPa•mm.The length parameter is 3 times larger than the characteristic element size, l = 0.06 mm.The discretization is enhanced in regions where crack propagation is anticipated (Fig. 27b), resulting in a mesh comprising 50,000 elements with a characteristic element length of h=0.02 mm along the crack path.Computations are conducted in a displacement-driven context, with the displacement increment adjusted to u = 110 -3 mm as the peak load is approached.The resulting load-displacement curve according to the 3D phase-field model is given in Fig. 27c.Comparing 3D simulation results (Fig. 19c) with simple 3D plane problem with through-thickness crack, it can be noted that the predicted profile of the loading-displacement relation for two surface flaws show similar behavior, nevertheless, crack interacting process (Fig. 27c) yields a moderate smooth softening behavior before the final failure.The crack propagation contours at various stages of loading in terms of phase fields for the two preexisting semi-circular surface flaws are given in Fig. 28.Four stages of the cracking process are shown, from the initiation of damage to the complete rupture of the plate.Blue and red colors correspond to the completely intact and the fully broken state of the material, respectively.First, damage initiates at adjacent nearby corner points of the front of each of the two cracks and propagates toward each other (Fig. 28a).This stage of fracture demonstrates the interaction of two separate preexisting semicircular surface cracks.Then, cracking initiates along the entire perimeter of each front and eventually both propagating cracks coalescence with each other combining into one common defect (Fig. 28b).After coalescence at the center of the plate and already the general semi-elliptical crack continues to propagate (Fig. 28c).According to Fig. 28c, the crack initiation takes place after reaching the peak load at values of applied displacements u = 0.0098mm while the complete failure occurs with the formation of a through-thickness crack (Fig. 28d) at displacements u = 0.014mm.These results demonstrate that the phase field model can inherently predict crack propagation caused by surface damage under monotonic loading conditions in three-dimensional contexts, regardless of the geometry and dimensions involved.

Figure 2 :
Figure 2: Phase field fracture scanning electron microscope images: (a) general field, (b) main crack tip detail.

Figure 4 :Figure 5 :Table 3 :
Figure 4: Single-edge-notched tension test: (a) boundary conditions and geometry, (b) finite element mesh.Representative results for notched square plate subjected to tension test shows in Fig.5(a) the force versus displacement curve.As observed in the figure, damage brings in an important drop in the load, with the crack propagating along the symmetry plane of the specimen.The resulting fracture patterns at different steps are illustrated in Figs.5(b-d).Blue and red colors correspond to the completely intact and the fully broken state of the material, respectively.A quantitative agreement with the results of Miehe et al.[3] and Martínez-Pañeda and Golahmar[9] was observed for the same governing parameters of computations.

Figure 6 :Figure 7 :
Figure 6: Force versus displacement curve (a) and ductile fracture pattern of a cracked plate at a displacement of (b) u = 0.002 mm, (c) u = 0.01 mm.Similar behavior of phase fields in the simple benchmark of a cracked square plate subjected to tension has been discussed in the literature.According to Martínez-Pañeda et al.[9], if the phase field fracture driving force is taken to be the elastic and plastic part of the total strain energy density, then cracking patterns can follow plastic localisation regions like a slip band, depending on the boundary value problem and the material properties (yield stress, G c , etc.).Obviously, ductile phase field fracture pattern will be observed in pure mode I if the configuration and loading conditions of a compact-tension sample are simulated.Next we considered the finite element computations of the crack phase-field  in the notched square plate (Fig.4) subjected to tension test according to the simple benchmark procedure for given values of the of the material fracture toughness G c length-scale parameter l and mesh size h.Specifically, Fig.7(a) shows the force versus displacement curves obtained for three different values of the critical energy release rate G c with fixed values of other parameters (l =0.04, h=0.01, u = 110 - 5 mm) for brittle fracture.As expected given the choice of G c = 1.3, 2.7 and 5.4 MPa•mm leads to increase in the values of maximum loads at fracture.Fig.7(b) shows the force versus displacement curves for ductile fracture.In this case, an increase of the G c value leads to an increase of the critical strains.

Figure 11 :
Figure 11: Single-edge-notched shear test.Crack patterns for  =0˚ are depicted on the left-hand side at the deformation state: (a) u =0.008 mm; (b) u =0.0125.Results obtained for  =7˚ are given on the right-hand side at the deformation state: (c) u =0.007;(d) u=0.0115 mm.

Figure 14 :
Figure 14: Force versus displacement curves (a) and ductile fracture pattern at a displacement of (b) u =0.138 mm, (c) u=0.300 mm.

Figure 17 :Figure 18 :
Figure 17: Experimental crack paths (a,c,e) and predicted fracture patterns (b,d,f) for pure mode I (a,b) and initial pure shear (c-f).

Figure 22 :
Figure 22: Configuration of biaxially loaded plate with a semi-elliptical surface crack: (a) representative volume; (b) global finite element mesh.

Figure 23 :
Figure 23: Load-displacement curves as a function of loading biaxiality for aspect ratio (a) a/c = 1 (semi-circle) and (b) a/c = 0.5 (semiellipse).The force versus displacement responses as a function of nominal stress biaxial ratio  =  xx / yy for aspect ratio a/c=1 (semi-circle) and a/c=0.5 (semi-ellipse) are shown in Fig.23.As stated in Shlyannikov et al. the loading biaxiality has a strong influence on the material property[27].As anticipated, a plate initially containing a surface defect exhibits greater loadbearing capacity under uniaxial tension compared to biaxial loading conditions.Specifically, for plates with pre-existing surface defects, the stiffness of the material degrades more rapidly when the defect is semi-elliptical rather than semi-circular,

Figure 24 :Figure 25 :Figure 26 :
Figure 24: Contour plots of the phase-field as the surface crack progresses through mutually perpendicular sections for pure mode I fracture under equi-biaxial tension: (a) u = 0.0029 mm, (b) u = 0.0063 mm, (c) u = 0.0069 mm.

Figure 27 :
Figure 27: Three-dimensional mode-I tension test with two surface cracks: (a) geometry, loading and boundary conditions, (b) finite element mesh, (c) load-displacement curve.

Table 1 :
ANSYS DOF command labels and corresponding forcing quantities.

Table 2 :
Main mechanical properties and loading conditions.

Table 4 :
Main elastic-plastic mechanical properties of the material and loading conditions.

Table 5 :
Main elastic-plastic mechanical properties of the material and loading conditions.

Table 6 :
Main mechanical properties and loading conditions.

Table 7 :
Main mechanical properties and loading conditions.