# An open-source implementation of a phase-field model for brittle fracture using Gridap in Julia

This article proposes an open-source implementation of a phase-field model for brittle fracture using a recently developed finite element toolbox, Gridap in Julia. The present work exploits the advantages of both the phase-field model and Gridap toolbox for simulating fracture in brittle materials. On one hand, the use of the phase-field model, which is a continuum approach and uses a diffuse representation of sharp cracks, enables the proposed implementation to overcome such well-known drawbacks of the discrete approach for predicting complex crack paths as the need for re-meshing, enrichment of finite element shape functions and an explicit tracking of the crack surfaces. On the other hand, the use of Gridap makes the proposed implementation very compact and user-friendly that requires low memory usage, and provides a high degree of flexibility to the users in defining weak forms of partial differential equations. A test on a notched beam under symmetric three-point bending and a set of tests on a notched beam with three holes under asymmetric three-point bending is considered to demonstrate how the proposed Gridap based phase-field Julia code can be used to simulate fracture in brittle materials.

There are no comments yet.

## Authors

• 1 publication
• ### Fireshape: a shape optimization toolbox for Firedrake

We introduce Fireshape, an open-source and automated shape optimization ...
05/14/2020 ∙ by Alberto Paganini, et al. ∙ 0

• ### OptiMic: A tool to generate optimized polycrystalline microstructures for materials simulations

Polycrystal microstructures, with their distinct physical, chemical, str...
07/06/2021 ∙ by Prince Henry Serrao, et al. ∙ 0

• ### A simple and robust Abaqus implementation of the phase field fracture method

The phase field fracture method is attracting significant interest. Phas...
04/16/2021 ∙ by Yousef Navidtehrani, et al. ∙ 0

• ### Tusas: A fully implicit parallel approach for coupled nonlinear equations

We develop a fully-coupled, fully-implicit approach for phase-field mode...
06/27/2020 ∙ by Supriyo Ghosh, et al. ∙ 0

• ### Multidimensional coupling: A variationally consistent approach to fiber-reinforced materials

A novel mathematical model for fiber-reinforced materials is proposed. I...
01/07/2021 ∙ by Ustim Khristenko, et al. ∙ 0

• ### T-IFISS: a toolbox for adaptive FEM computation

T-IFISS is a finite element software package for studying finite element...
08/15/2019 ∙ by Alex Bespalov, et al. ∙ 0

• ### A unified Abaqus implementation of the phase field fracture method using only a user material subroutine

We present a simple and robust implementation of the phase field fractur...
04/08/2021 ∙ by Yousef Navidtehrani, et al. ∙ 0

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

To design structures with high reliability, failure analysis of structures is of great importance in engineering applications. As fracture due to crack initiation and propagation is one of the most often encountered failure modes in engineering materials and structures, modeling of fracture in solids has always been one of the most intriguing topics of research interests. Numerical modeling of fracture in solids has mainly been done either by using a discrete or a continuum approach. In the discrete approach, cracks in the material body are modeled as the discontinuity of the displacements in the domain whereas in the continuum approach a diffused approximation of cracks is used to model fracture as a continuum damage process for which displacements are continuous but the material stiffness gradually degrades. Linear elastic fracture mechanics (LEFM) griffitli1920phenomena; irwin1956onset; williams2001introduction; luo2018linear and cohesive zone model (CZM) barenblatt1962mathematical; dugdale1960yielding are the notable theories in the category of the discrete approach for fracture modelling. Although LEFM and CZM are very popular, their implementation requires an explicit tracking of the discontinuity in the displacement field that poses difficulty in modeling an arbitrary complex crack path.

Knowing the well-known drawbacks of the discrete approach for modeling complicated crack paths, researchers generally refer to the continuum approach that provides the crack paths as part of the solutions of the governing partial differential equations. One of the most popular theories in the category of continuum approach is the phase-field model (PFM) ambati2015review. There are of course several phase-field approaches to brittle fracture that have been independently developed in the mechanics community francfort1998revisiting; bourdin2000numerical; bourdin2008variational; kuhn2008phase; amor2009regularized; kuhn2010continuum; miehe2010phase; borden2014higher; dhas2018phase and in the physics community aranson2000continuum; karma2001phase; hakim2009laws; spatschek2011phase; eastgate2002fracture; henry2004dynamic as well. In this article, a phase-field model proposed by Dhas et al. dhas2018phase is adopted as the model provides a thermodynamically consistent way of accommodating dissipative energy effects whenever needed. In all the phase-field models, a diffused approximation of sharp cracks is used by introducing a length scale parameter and an internal variable called phase-field. The accuracy of the diffused approximation depends on the value of the length scale parameter and may represent the original crack problem if the length scale parameter is chosen sufficiently small. Although this feature of phase-field models imposes a highly efficient implementation of the model while using the finite element method as very fine meshes are required for regularizing the sharp cracks with a small value of length scale parameter, the model has gained huge popularity in the research community as it can be incorporated in commercial finite element software such as Abaqus msekh2015abaqus; liu2016abaqus; molnar20172d; wu2020comprehensive; navidtehrani2021simple; navidtehrani2021unified. However, to make the phase-field model available for a wider class of practitioners and researchers, there are also attempts towards open-source implementation of phase-field models by using finite element method natarajan2019fenics

and machine learning-based approaches

samaniego2020energy; goswami2020transfer.

In this article, a new open-source implementation of a phase-field model is proposed using a recently developed finite element toolbox Gridap available in the programming language Julia bezanson2012julia; bezanson2017julia that shares the advantages of both the static and dynamic languages. The programming language Julia is computationally efficient as the static languages such as Fortran, C++, etc., and also easy to use as the dynamic languages like Matlab, Python, Mathematica, etc. Gridap is an extensible finite element toolbox badia2020gridap; verdugo2019user in Julia that can be used to solve a wide range of physical problems modeled mathematically using partial differential equations (PDEs). In contrast to other finite element libraries written in Julia such as FinEtools, JuAFEM, and JuliaFEM frondelius2017juliafem, Gridap uses a novel software design (for example, high-level API calls) that enables one to compute the value for a specific cell on the fly and never store the values for all cells in the mesh simultaneously and thus essentially requires very low memory usage. Moreover, Gridap provides a high degree of flexibility to the users as they can implement any PDEs-based mathematical model such as a phase-field model using a very compact syntax without explicitly writing any for-loop for assembly over elemental matrices. To develop a Gridap based open source program for phase-field modeling of brittle fracture, a thermodynamically consistent phase-field model is briefly described first in Section 2. Then, the derivation of the weak form corresponding to the governing PDEs of the phase-field model and the finite element implementation in Julia are provided in Section 3. Successful implementation of the phase-field model using Gridap is demonstrated in Section 4 through a test on a notched beam under symmetric three-point bending and a set of tests on a notched beam with three holes under asymmetric three-point bending tests. Finally, the outcomes of the present work are summarized, and concluding remarks are accordingly made in Section 5.

## 2 Phase-field model

In this section, a brief description of a thermodynamically consistent phase-field approach dhas2018phase to brittle fracture in elastic solids under small strains and isothermal conditions is provided.

### 2.1 Kinematics

Consider an open set to be the reference configuration of a deformable body in the three dimensional Euclidean space . Let and be the smooth boundary and the closure of , respectively. The displacement field may be defined as at any instant of time

. Within a small deformation set-up, the strain tensor

is given by

 ϵ(u)=12(∇u+∇uT), (1)

where denotes transpose of a tensor and the gradient operator. In the phase-field model, sharp cracks are approximated by using a diffused representation of the cracks via a length scale parameter and an internal variable called phase-field (see Fig. 1). Using the values of the phase-field variable, one can describe the damaged, undamaged or partially damaged states of matter as follows: for the undamaged state, for the fully damaged state and for a partially damaged state. Considering fully damaged states as the fracture, crack set may be defined as . Here, the damage process is considered to be irreversible that is if at time then for all . Thus, the deformed and damaged states of the material body may be described by considering phase-field variable

as an additional kinematic descriptor along with the displacement vector

.

### 2.2 Force Balances

To describe the deformation of a material body under external loading, one can derive the force balances from a virtual power principle by considering a macro- and a micro-force system gurtin1996generalized. While stress tensor , traction vector and body force define the macro-force system, the micro-force system includes scalar micro-traction , vector micro-stress and scalar micro-stress . It is important to note that the micro-system is also at the same continuum level of the macro-system and the coinage ‘micro-force’ may be a misnomer. One can identify as the power conjugate of and , as the power conjugate of and express the external power for any arbitrary part of the body as

 Pext=∫∂P[t(n)⋅˙u+χ(n)⋅˙s]dA+∫Pb⋅˙udV, (2)

where denotes the time derivative of a variable, the unit normal vector to the boundary , and are, respectively, measures on and . The internal power can be defined by the summation of power expenditure of over , over and over , and can be given by

 Pint=∫P(σ:∇˙u+ξ⋅∇˙s+π˙s)dV. (3)

Denoting the virtual counterparts of and by and , respectively, one can define a set called as the generalized virtual velocity vector. Then the external virtual power and the internal virtual power can be expressed as

 ~Pext(P;V)=∫∂P[t(n)⋅~u+χ(n)⋅~s]dA+∫Pb⋅~udV, (4)

and

 ~Pint(P;V)=∫P(σ:∇~u+ξ⋅∇~s+π~s)dV, (5)

respectively. Employing Eq. (4), Eq. (5) and invoking the virtual power principle i.e., , one can arrived at

 ∫∂P[t(n)⋅~u+χ(n)⋅~s]dA+∫Pb⋅~udV=∫P(σ:∇~u+ξ⋅∇~s+π~s)dV. (6)

Appropriate choices of in Eq. (6) may lead to the macro- and the micro-force balances as described below.

#### 2.2.1 Macro-force balance

Considering i.e., by substituting in Eq. (6), the macro-force balance equation may be obtained as

 ∫∂P(t(n)−σn).~udA=−∫P(∇⋅σ+b)⋅~udV, (7)

which holds for all and any arbitrary sub-domain . Applying the localization theorem on Eq. (7), one can get that

 t(n)=σn, (8)

which is the macro-traction condition and

 ∇⋅σ+b=0. (9)

called the macro-force balance. Equations (8) and (9) may be identified as the classical traction condition and the local linear momentum balance equation, respectively.

#### 2.2.2 Micro-force balance

Considering i.e., by substituting in Eq.(6), the equation for micro-force balance can be derived as

 ∫∂P[χ(n)−ξ⋅n]⋅~sdA=−∫P(∇⋅ξ−π)⋅~sdV. (10)

Since Eq. (10) holds for all and any sub-domain , Eq. (10) may be localized as

 χ(n)=ξ⋅n, (11)

and

 ∇⋅ξ−π=0. (12)

Equations (11) and (12) are called the micro-traction condition and the micro-force balance, respectively.

#### 2.2.3 Thermodynamics and constitutive modeling

In this section, constitutive relations for the macro- and micro- stresses are derived by imposing the first and second laws of thermodynamics. Considering an iso-thermal condition, one may state the second law of thermodynamics for any sub-domain as the free energy inequality:

 ddt∫PψdV≤Pext, (13)

where is the Helmholtz free-energy of the system. Using Eq. (2), Eq. (3) and the power balance, i.e. = , the inequality given by Eq. (13) may be expressed as

 ∫P˙ψdV≤∫P(σ:˙ϵ+ξ⋅∇˙s+π˙s)dV, (14)

where the equality (which follows from the symmetry of ) is used. Since, the inequality given by Eq. (14) holds for any arbitrary sub-domain , one can have

 ˙ψ−(σ:˙ϵ+ξ⋅∇˙s+π˙s)≤0, (15)

which must be satisfied whilst determining or postulating the constitutive relations for the thermodynamic fluxes , and in terms of the kinematic quantities , and .

### 2.3 Constitutive response functions

Let the free energy of the system be function of , , and may be written as

 ψ=^ψ(ϵ,∇s,s). (16)

Using the chain rule in Eq. (

16), one can get the rate of free energy as

 ˙ψ=∂ϵψ:˙ϵ+∂∇sψ⋅∇˙s+∂sψ˙s, (17)

where with a suffix represents the derivative of a function with respect to the argument in the suffix while keeping others fixed. Considering that the scalar micro-stress has an energetic part and a dissipative part i.e., and substituting Eq. (17) in Eq. (15), one can get that

 (σ−∂ϵψ):˙ϵ+(ξ−∂∇sψ)⋅∇˙s+(πen−∂sψ)˙s+πdis˙s≥0. (18)

Applying the Coleman-Noll procedure coleman1974thermodynamics to Eq. (18), one can arrive at the constitutive relations for the thermodynamics fluxes as

 σ=∂ϵψ, (19)
 ξ=∂∇sψ, (20)
 πen=∂sψ, (21)

 πdis˙s≥0. (22)

Determination of the constitutive relation for must be done in such a way that the inequality constraint given by Eq. (22) always satisfy. From the irreversiblity condition on damage i.e. , it can be seen that a possible choice for could be , where is a constitutive function with .

### 2.4 Specialized constitutive relations

To quantify the thermodynamic forces, one need to specialize the constitutive relations by postulating an explicit expression of the Helmholtz free energy in terms of , and . Let the Helmholtz free energy be a sum of elastic energy and fracture energy as

 ψ(ϵ,∇s,s)=ψelas(ϵ,s)+ψfrac(∇s,s). (23)

It is assumed that crack cannot propagate under pure compression and imposed by considering an additive decomposition of into a volumetric part and a deviatoric part i.e.

 ϵ=ϵvol+ϵdev, (24)

where

 ϵvol=Pvolϵ;ϵdev=Pdevϵ. (25)

In Eq. (25), volumetric and deviatoric parts of a second order tensor are obtained by introducing fourth order projection tensors and , respectively. Defining with denoting the fourth order elasticity tensor, the elastic part of the free energy may be postulated as

 ψelas(ϵ)=s2ψelas+(ϵ)+ψelas−(ϵ), (26)

where is the elastic energy part due to a combination of pure tension and shear,

 ψelas+(ϵ)=12(⟨p⟩+I:ϵvol+PdevCϵ:ϵ% dev), (27)

and is the elastic energy part due to pure compression

 ψelas−(ϵ)=12⟨p⟩−I:ϵvol. (28)

In Eq. (26), a degradation function is introduced to account for the reduced stiffness of the material due to damage. Note that the degradation function is only associated with the so-called positive part of the elastic energy to impose the condition that cracks cannot propagate under pure compression. In Eq. (27), and in Eq. (28) . The fracture energy may be postulated as

 ψfrac(∇s,s)=Gc((1−s)22ls+ls2∇s⋅∇s), (29)

where is the critical energy release rate and is the phase-field length scale parameter. Employing equations (19), (20), (21), (27), (28) and (29), explicit expressions for the thermodynamic fluxes , and may be derived as

 σ=∂ϵψ=s2(⟨p⟩+I+PdevCϵ)+⟨p⟩−I, (30)
 ξ=∂∇sψ=Gcls∇s, (31)

and

 πen=∂sψ=2sψelas+(ϵelas)−Gcls(1−s), (32)

respectively. Substituting the above constitutive relations in the macro- and micro-force balances, one may express the governing partial differential equations in terms of the kinematic descriptors and .

### 2.5 Boundary value problem and the strong form

Using the expressions derived in the previous section, one may write the strong form of the governing PDEs as a boundary value problem and derive the corresponding weak form of the governing PDEs for the finite element formulation of phase-field model. One can express the stress tensor in terms of kinematic variables as by defining as

 Cmod={s2C,forp≥0s2(PdevC)+PvolC,forp<0. (33)

Using , Eq. (9) may be re-written as

 ∇⋅(Cmodϵ)+b=0. (34)

Similarly, using the expressions of and and assuming , Eq. (12) may be re-written as

 ∇⋅(Gcls∇s)−2sψelas+(ϵ)+Gcls(1−s)=0. (35)

To account for the irreversiblity condition on damage, a history function , where for any input argument , of the so-called positive part of elastic energy is employed miehe2010phase. Using the history function , Eq. (35) may be expressed as

 ∇⋅(Gcls∇s)−2sH(E)+Gcls(1−s)=0. (36)

The governing PDEs (34) and (36) are coupled and subject to boundary conditions, such as prescribed displacement and applied traction on and , respectively. Equations (34) and (36) together with the boundary conditions are called the strong form of the governing PDEs. In absence of body force i.e. , the strong form for phase-field modeling of brittle fracture in an elastic solid defined by domain subjected to displacement boundary condition on the boundary can be given in a compact form as

 ∇⋅σ=0(% whereσ=Cmodϵ)inΩ, (37a) u=¯uon∂Ωu, (37b) σn=0on∂Ω∖∂Ωu, (37c)

corresponding to the macro-system and

 ∇⋅(GclA∇s)−2sH(E)+Gcl(1−s)=0inΩ, (38a) ∇s⋅n=0on∂Ω∖∂Ωu, (38b)

corresponding to the micro-system. In the present study, beams are made of isotropic materials for which components of the fourth order elasticity tensor may be expressed as , where and are the Lamé parameters and denotes the Kronecker delta. The Lamé parameters are related to Young’s modulus and Poisson’s ratio by and .

## 3 The weak form and an outline for finite element implementation in Julia

In this section, first the weak form for the strong form given by Eq. (37) and Eq. (38) is derived and then an outline for the finite element formulation is provided through a numerical example. Let and be the test functions corresponding to displacement and phase-field , respectively. The trial spaces with the given displacement boundary condition may be given by

 Hu={u∈H1(Ω);u=¯uon∂Ωu}, (39)
 Hs={s∈H1(Ω)}, (40)

where is a prescribed displacement on . The test spaces may be defined as

 Vu={v∈H1(Ω);v=0on∂Ωu}, (41)
 Vs={ϕ∈H1(Ω)}. (42)

One can obtain the following weak form: Find and such that for all and ,

 aDisp(u,v)=bDisp(v), (43a) aPF(s,ϕ)=bPF(ϕ), (43b)

where

 aDisp(u,v)=∫Ωϵ(v):σ(u)dΩ,bDisp(v)=0, (44a) aPF(s,ϕ)=∫Ω(Gcls∇s⋅∇ϕ+2sϕH(E)+Gclssϕ)dΩ,bPF(ϕ)=∫ΩGclsϕdΩ. (44b)

In the present study, a staggered scheme originally proposed by Miehe et al. miehe2010phase is employed to solve for the unknown displacement vector and phase-field from the weak form defined by Eq. (43) and Eq. (44) using Gridap.

One of the salient features of Gridap is that one can directly use the weak form in Julia for the finite element implementation using Gridap. All the steps associated with the finite element simulations such as creating the mesh file, implementing the weak form, application of boundary conditions, solutions for the unknown field variables, and post-processing for the output files are described through a Julia code on numerical simulation of a test on a notched beam under symmetric three-point bending (see Section 4.1 for numerical results). One may readily apply the developed phase-field-based Julia code for simulating other brittle fracture problems with appropriate modifications. For reproducing the results presented in Section 4, one needs to first load the following Julia packages in the script file which is presently written in a jupyter notebook.

using GridapGmsh
using Gridap
using Gridap.Geometry
using Gridap.TensorValues
using PyPlot

One can define the input parameters associated with the elastic and fracture material properties of the notched beam by writing the following lines.

const E_mat = 20.8
const _mat = 0.3
const Gc = 5e-4
const ls = 0.03
const  = 1e-15

For the finite element simulations, one needs to have the discretization of the computational domain which can be generated by writing a mesh file in Julia (see A) and load that mesh file to build an instance of “DiscreteModel” by the following lines.

model = GmshDiscreteModel("BeamWithNotchSymThreePtBending.msh")
writevtk(model,"BeamWithNotchSymThreePtBending")

For an isotropic material, the constitutive tensor can be defined for two-dimensional plane stress and plane strain problems by writing the following function.

function ElasFourthOrderConstTensor(,,PlanarState)
# 1 for Plane Stress and 2 Plane Strain Condition
if PlanarState == 1
C1111 = /(1-*)
C1122 = (*)/(1-*)
C1112 = 0.0
C2222 = /(1-*)
C2212 = 0.0
C1212 = /(2*(1+))
elseif PlanarState == 2
C1111 = (*(1-*))/((1+)*(1--2**))
C1122 = (*)/(1--2**)
C1112 = 0.0
C2222 = (*(1-))/(1--2**)
C2212 = 0.0
C1212 = /(2*(1+))
end
C_ten = SymFourthOrderTensorValue(C1111,C1112,C1122,C1112,C1212,C2212,C1122,C2212,C2222)
return  C_ten
end

In the present study, plane strain condition is assumed for which the constitutive tensor is computed by calling the above Julia function as given below.

const C_mat = ElasFourthOrderConstTensor(E_mat,_mat,2)

To satisfy the assumption that crack can not propagate under pure compression, stress and strain tensors are decomposed into a volumetric and a deviatoric part by introducing the projection operators and , respectively, which are defined by the following lines.

I2 = SymTensorValue{2,Float64}(1.0,0.0,1.0)
I4 = I2$\otimes$I2
I4_sym = one(SymFourthOrderTensorValue{2,Float64})
P_vol = (1.0/3)*I4
P_dev = I4_sym - P_vol

To express the stress tensor in terms of kinematics variables i.e. , where the expression of is given by Eq. (33), the following function is defined in Julia.

function fun(,_in,s_in)
_elas = C_mat$\odot$$\varepsilon if tr(_in) >= 0 = (s_in ^2+)*_elas elseif tr(_in) < 0 = (s_in ^2+)*P_dev _elas + P_vol _elas end return end One can determine the so-called positive part of elastic free energy, which is given by Eq. (27), by writing the following function in Julia. function Pos(_in) _elas = C_mat\odot$$\varepsilon$_in
if tr(_in) >= 0
Plus = 0.5*(_in  _elas)
elseif  tr(_in) < 0
Plus = 0.5*((P_dev  _elas)(P_dev$\odot$$\varepsilon$_in))
end
return Plus
end

One needs to generate a discrete approximation of the finite element test and trial spaces of the problem on the discretized computational domain. Approximation of the finite element spaces associated with the phase field variable can be done by the following lines in Julia.

order = 1
reffe_PF = ReferenceFE(lagrangian,Float64,order)
V0_PF = TestFESpace(model,reffe_PF;
conformity=:H1)
U_PF = TrialFESpace(V0_PF)
sh = zero(V0_PF)

Similarly, one can generate the approximation of the finite element spaces associated with the displacement variable by writing the following lines in Julia.

reffe_Disp = ReferenceFE(lagrangian,VectorValue{2,Float64},order)
V0_Disp = TestFESpace(model,reffe_Disp;
conformity=:H1,
uh = zero(V0_Disp)

To compute the integrals in the weak form given by Eq. (44) numerically, one needs to define an integration mesh along with an integration rule (for example, Gauss quadrature) in each of the cells in the triangulation. Using Gridap, one can easily define the integration mesh and the corresponding Lebesgue measure by using the built-in functions “Triangulation” and “Measure”, respectively. For instance, one can use the following lines for integrating the weak form given by Eq. (44) defined on the domain

using a quadrature rule of degree two times the order of interpolation in the cells of the triangulation.

degree = 2*order
= Triangulation(model)
d$\Omega$ = Measure(,degree)

One can determine the applied load on a part of the boundary of the domain by determining boundary integral using the following built-in functions available in Gridap.

labels = get_face_labeling(model)
d$\Gamma$_Load = Measure(_Load,degree)
n_$\Gamma$_Load = get_normal_vector(_Load)

To find the values of variables that are defined by using the built-in function “CellState” at the Gauss points, one need to use the “project” function as defined below.

function project(,model,d$\Omega$,order)
reffe = ReferenceFE(lagrangian,Float64,order)
V = FESpace(model,reffe,conformity=:L2)
a(,) = (*)*d$\Omega$
b() = (*)*d$\Omega$
op = AffineFEOperator(a,b,V,V)
qh = solve(op)
return qh
end

In the present study, a staggered scheme is used to update the solution from the pseudo time to miehe2010phase. Given the displacement vector, phase-field and the history function at the time , one can update the phase-field at the time by using the following function.

function  stepPhaseField(uh_in,PlusPrev_in)
a_PF(,) = ( Gc*ls* + 2*PlusPrev_in** + (Gc/ls)** )*d$\Omega$
b_PF() = ( (Gc/ls)*)*d$\Omega$
op_PF = AffineFEOperator(a_PF,b_PF,U_PF,V0_PF)
sh_out = solve(op_PF)
return sh_out
end

Using the values of displacement vector and the history function at time , and the computed value of phase-field at time , one can update the displacement vector at time by calling a function as given below.

function   stepDisp(uh_in,sh_in,vApp)
uApp1() = VectorValue(0.0,0.0)
uApp2() = VectorValue(0.0,0.0)
uApp3() = VectorValue(0.0,-vApp)
U_Disp = TrialFESpace(V0_Disp,[uApp1,uApp2,uApp3])
a_Disp(,) = ( ( (fun$\circ$(,(uh_in),sh_in)) ) )*d$\Omega$
b_Disp() = 0.0
op_Disp = AffineFEOperator(a_Disp,b_Disp,U_Disp,V0_Disp)
uh_out = solve(op_Disp)
return uh_out
end

Once both displacement and phase field are determined at time , one can update the energy history function at time by defining it as

 Hn+1={En+1forEn+1>EnEnotherwise, (45)

where the history function of elastic free energy need to be defined in Julia (see Section 2.5), which can be achieved by writing the following lines.

function new_EnergyState(PlusPrev_in,hPos_in)
Plus_in = hPos_in
if Plus_in >= PlusPrev_in
Plus_out = Plus_in
else
Plus_out = PlusPrev_in
end
true,Plus_out
end

Finally, one can simulate the brittle fracture in the notched beam by applying a monotonic displacement control loading and solve for the unknown displacement and phase-field at each loading step by writing the main routine in Julia that uses the above-defined functions as listed below.

vApp = 0
delv = 1e-3
const vAppMax = 0.1
innerMax = 10
count = 0
Displacement = Float64[]
push!(Displacement, 0.0)
sPrev = CellState(1.0,d$\Omega$)
sh = project(sPrev,model,d$\Omega$,order)
PlusPrev = CellState(0.0,d$\Omega$)
while vApp .< vAppMax
count = count .+ 1
if vApp >= 3e-2
delv = 1e-4
end
vApp = vApp .+ delv
print("\n Entering  displacemtent  step :", float(vApp))
for inner = 1:innerMax
hPlusPrev = project(PlusPrev,model,d$\Omega$,order)
RelErr = abs(sum(( Gc*ls*(sh)(sh) + 2*hPlusPrev*sh*sh + (Gc/ls)*sh*sh)*d$\Omega$-( (Gc/ls)*sh)*d$\Omega$))/abs(sum(( (Gc/ls)*sh)*d$\Omega$))
sh = stepPhaseField(uh,hPlusPrev)
uh = stepDisp(uh,sh,vApp)
hPos_in = Pos$\circ$((uh))
update_state!(new_EnergyState,PlusPrev,hPos_in)
if  RelErr < 1e-8
break
end
end
Node_Force = sum((n_$\Gamma$_Load$\cdot$(fun$\circ$((uh),(uh),sh))) *d$\Gamma$_Load)
push!(Displacement, vApp)
end
end

One can create output files at each loading step that can be viewed in ParaView. For instance, one can create a “.vtu” file to save data for the solution of the displacement vector and phase-field at each loading step by using the following lines in Julia.

writevtk(,"results_SymThreePtBendingTest",cellfields
= ["uh"=>uh,"s"=>sh])

One can generate the load-displacement curve by using the plot command as given below.

## 4 Numerical Simulations

Proposed phase-field-based Julia codes are validated against a test on a notched beam under symmetric three-point bending and a set of tests on a notched beam with three holes under asymmetric three-point bending. The effect of length scale parameter value on the fracture response is well understood and hence not repeated here. To validate the proposed open-source implementation, one particular value of the length scale parameter , which is typically mentioned in the literature for the given problem, is taken. For the numerical simulation, a non-uniform finite element mesh with a finer mesh (length of the largest side of triangular elements is less than half of value) in regions where cracks may propagate is used. For finite element mesh generation of the notched beams used for symmetric and the asymmetric three-point bending test, Julia codes are provided in A and B, respectively.

### 4.1 Symmetric three point bending test

Modeling of brittle fracture in a simply supported notched beam under symmetric three-point bending is one of the classical benchmark problems which has frequently been analyzed in the literature ambati2015review; miehe2010phase; miehe2007robust; wu2018length. The three-point bending test set-up and a finite element mesh used for the simulation are demonstrated in Fig. 2.

For the numerical simulation, material properties for the notched beam are taken as , , and . Displacement control loading (monotonic displacement is applied in small increments ) is considered and the damage profiles for the notched beam at different stages of applied displacement are presented in Fig. 3.

For applied displacement up to , monotonic increment of and from to until failure (here ) monotonic increment of is used. As can be seen from Fig. 4, load displacement curve using the proposed open source implementation matches quite well with the results reported in the literature miehe2010phase.

### 4.2 Asymmetric notched three point bending test

The developed Julia code for the phase-field model is validated against a set of tests on a notched beam with three holes under asymmetric three-point bending, which were carried out by Ingraffea and Grigoriu ingraffea1990probabilistic and numerically analyzed in Bittencourt et al. bittencourt1996quasi. Material parameters are taken as , , and . To verify whether the proposed Julia implementation of phase-field model can predict experimentally observed complex crack patterns, three different conﬁgurations of the specimen characterized by the values and (see Fig. 5 for geometry and the boundary conditions, and a finite element mesh used for simulation) are considered.

Prediction of crack path for the beam with three holes and a pre-notch defined by (a) and , (b) and and (c) and are considered. Displacement control loading (monotonic displacement is applied in small increments ) is considered and the damage profiles for the beam with three holes and a pre-notch defined by (a) and , (b) and and (c) and at different stages of applied displacement are presented in Fig. 6, Fig. 8 and Fig. 10, respectively. As can be seen from Fig. 7, Fig. 9 and Fig. 11, the proposed Julia implementation shows a very good prediction of the experimentally observed crack paths which are very sensitive to the height and relative location of the pre-notch. Remarkably, the proposed implementation reproduces the intricate deviation of crack path due to the local stress concentration around the bottom hole which is experimentally observed (see Fig. 9).

## 5 Concluding Remarks

The present study has provided a novel numerical implementation of a thermodynamically consistent phase-field model for brittle fracture using an open-source ﬁnite element toolbox, Gridap in Julia. The proposed implementation is validated against a few numerical and experimental results available in the literature. As the proposed implementation is available with an open-source license, it may eliminate the technical barrier for practitioners and researchers who are interested to explore the phase-field model for solving a wide range of brittle fracture problems. Moreover, the proposed implementation will expose the users to many built-in packages of Julia that may be useful for researchers who want to extend the proposed implementation for the case of ductile fracture or other applications. Most importantly, the availability of an open-source code that is compact, user friendly, highly efficient, and accessible to everyone will allow a third-party verification and essentially establish a high standard for efficient open-source code development.

## Data accessibility

The present work does not generate any experimental data. Julia scripts used as the source codes for symmetric three-point bending test are provided in Section 3 of this article. Jupyter notebook files for the proposed Julia implementation of a phase-field model can be downloaded from Julia code for phase-field model.

## Declaration of Competing Interest

The author of this article declares that he has no conflict of interest.

## Acknowledgments

The author of this article gratefully acknowledge support from the Indian Institute of Technology Bhubaneswar under the grant SP107.

## Appendix A Julia code to create finite element mesh file for a notched beam under symmetric three point bending

One can generate finite element mesh file in Julia by using the GMSH mesh generator, which can be loaded in Julia by writing the following line.

using Gmsh: gmsh

One can create the mesh file “BeamWithNotchSymThreePtBending.msh” using the following lines in Julia.

const L = 8.0
const LL = 0.475.*L
const LR = 0.525.*L
const H = 2.0
const CH = 0.4 #Crack height
const CW = 0.2 #Crack Width
const ls = 0.03
const hfc = ls/100 #Mesh size parameter
const hf = ls/2.1 #Mesh size parameter
const h = 100*hf #Mesh size parameter
gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 1)
gmsh.model.geo.addPoint((L/2)+(CW/2), 0.0 , 0.0, hf,1)
gmsh.model.geo.addPoint(L, 0.0, 0.0, h, 2)
gmsh.model.geo.addPoint(L, H, 0.0, h, 3)
gmsh.model.geo.addPoint(LR, H, 0.0, hf, 4)
gmsh.model.geo.addPoint(LL, H, 0.0, hf, 5)
gmsh.model.geo.addPoint(0.0, H, 0.0, h, 6)
gmsh.model.geo.addPoint(0.0, 0.0, 0.0, h, 7)
gmsh.model.geo.addPoint((L/2)-(CW/2), 0.0 , 0.0, hf,8)
gmsh.model.geo.addPoint((L/2), CH , 0.0, hfc, 9)
gmsh.model.setPhysicalName(2, 1, "Domain")
gmsh.model.setPhysicalName(0, 2, "LeftSupport")
gmsh.model.setPhysicalName(0, 3, "RightSupport")
gmsh.model.mesh.field.setNumber(10, "VIn", hf)
gmsh.model.mesh.field.setNumber(10, "VOut", h)
gmsh.model.mesh.field.setNumber(10, "XMin", (L/2)-CW)
gmsh.model.mesh.field.setNumber(10, "XMax", (L/2)+CW)
gmsh.model.mesh.field.setNumber(10, "YMin", 0)
gmsh.model.mesh.field.setNumber(10, "YMax", H)
gmsh.model.mesh.field.setAsBackgroundMesh(10)
gmsh.model.geo.synchronize()
gmsh.model.mesh.generate(2)
gmsh.write("BeamWithNotchSymThreePtBending.msh")
gmsh.finalize()

## Appendix B Julia code to create finite element mesh file for a notched beam with three holes under asymmetric three point bending

One can create the mesh file “AsymThreePtBending.msh” by using the following lines in Julia.

using Gmsh: gmsh
const L = 20.0
const LL = 0.475.*L
const LR = 0.525.*L
const H = 8.0
const CH = 1.5 #Crack height
const CW = L/2000 #Crack Width
const  = 5.15
const CP = L/2 -
const SD = 1.0
const HP = 6.0
const HR = 0.25
const HH1 = 2.75
const HH2 = 4.75
const HH3 = 6.75
const ls = 0.01
const hfc = ls/50  #Mesh size parameter
const hf = ls/2.1 #Mesh size parameter
const hfl = 50*hf  #Mesh size parameter
const hfh = hf  #Mesh size parameter
const h = 100*hf #Mesh size parameter
= /180
cr1 = CP+(CW/2) + HP*tan()
cr2 = CP-(CW/2) + HP*tan()
const FMR = 40*ls
gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 1)
p1 = gmsh.model.geo.addPoint(CP+(CW/2), 0.0 , 0.0, h)
p2 = gmsh.model.geo.addPoint(L-SD, 0.0, 0.0, h)
p3 = gmsh.model.geo.addPoint(L, 0.0, 0.0, h)
p4 = gmsh.model.geo.addPoint(L, H, 0.0, h)
p5 = gmsh.model.geo.addPoint(LR, H, 0.0, hfl)
p6 = gmsh.model.geo.addPoint(LL, H, 0.0, hfl)
p7 = gmsh.model.geo.addPoint(0.0, H, 0.0, h)
p8 = gmsh.model.geo.addPoint(0.0, 0.0, 0.0, h)
p9 = gmsh.model.geo.addPoint(SD, 0.0, 0.0, h)
p10 = gmsh.model.geo.addPoint(CP-(CW/2), 0.0, 0.0, h)
p11 = gmsh.model.geo.addPoint(CP-(CW/2), CH, 0.0, hfc)
p12 = gmsh.model.geo.addPoint(CP+(CW/2), CH, 0.0, hfc)
l1 = gmsh.model.geo.addLine(p1, p2)
l2 = gmsh.model.geo.addLine(p2, p3)
l3 = gmsh.model.geo.addLine(p3, p4)
l4 = gmsh.model.geo.addLine(p4, p5)
l5 = gmsh.model.geo.addLine(p5, p6)
l6 = gmsh.model.geo.addLine(p6, p7)
l7 = gmsh.model.geo.addLine(p7, p8)
l8 = gmsh.model.geo.addLine(p8, p9)
l9 = gmsh.model.geo.addLine(p9, p10)
l10 = gmsh.model.geo.addLine(p10, p11)
l11 = gmsh.model.geo.addLine(p11, p12)
l12 = gmsh.model.geo.addLine(p12, p1)
p13 = gmsh.model.geo.addPoint(HP-HR, HH1, 0.0, hfh)
p14 = gmsh.model.geo.addPoint(HP+HR, HH1, 0.0, hfh)
p15 = gmsh.model.geo.addPoint(HP, HH1, 0.0, hfh)
ca1 = gmsh.model.geo.addCircleArc(p13, p15, p14)
ca2 = gmsh.model.geo.addCircleArc(p14, p15, p13)
p16 = gmsh.model.geo.addPoint(HP-HR, HH2, 0.0, hfh)
p17 = gmsh.model.geo.addPoint(HP+HR, HH2, 0.0, hfh)
p18 = gmsh.model.geo.addPoint(HP, HH2, 0.0, hfh)
ca3 = gmsh.model.geo.addCircleArc(p16, p18, p17)
ca4 = gmsh.model.geo.addCircleArc(p17, p18, p16)
p19 = gmsh.model.geo.addPoint(HP-HR, HH3, 0.0, hfh)
p20 = gmsh.model.geo.addPoint(HP+HR, HH3, 0.0, hfh)
p21 = gmsh.model.geo.addPoint(HP, HH3, 0.0, hfh)
ca5 = gmsh.model.geo.addCircleArc(p19, p21, p20)
ca6 = gmsh.model.geo.addCircleArc(p20, p21, p19)
pg1 = gmsh.model.addPhysicalGroup(2, [ps1])
pg2 = gmsh.model.addPhysicalGroup(1, [l5])
pg3 = gmsh.model.addPhysicalGroup(0, [p9])
pg4 = gmsh.model.addPhysicalGroup(0, [p2])
gmsh.model.setPhysicalName(2, pg1, "Domain")
gmsh.model.setPhysicalName(0, pg3, "LeftSupport")
gmsh.model.setPhysicalName(0, pg4, "RightSupport")
p22 = gmsh.model.geo.addPoint(CP-(CW/2), 0.8*CH, 0.0, hf)
p23 = gmsh.model.geo.addPoint(CP+(CW/2), 0.8*CH, 0.0, hf)
p24 = gmsh.model.geo.addPoint(HP, cr1, 0.0, hf)
p25 = gmsh.model.geo.addPoint(HP, cr2, 0.0, hf)
l13 = gmsh.model.geo.addLine(p22, p24)
l14 = gmsh.model.geo.addLine(p23, p25)
gmsh.model.mesh.field.setNumbers(1, "EdgesList", [l13,l14])
gmsh.model.mesh.field.setNumber(2, "IField", 1)
gmsh.model.mesh.field.setNumber(2, "LcMin", hf)
gmsh.model.mesh.field.setNumber(2, "LcMax", h)
gmsh.model.mesh.field.setNumber(2, "DistMin", FMR)
gmsh.model.mesh.field.setNumber(2, "DistMax", 1.5*FMR)
l16 = gmsh.model.geo.addLine(p23, p5)
gmsh.model.mesh.field.setNumbers(3, "EdgesList", [l15,l16])
gmsh.model.mesh.field.setNumber(4, "IField", 3)
gmsh.model.mesh.field.setNumber(4, "LcMin", hfl)
gmsh.model.mesh.field.setNumber(4, "LcMax", h)
gmsh.model.mesh.field.setNumber(4, "DistMin", FMR)
gmsh.model.mesh.field.setNumber(4, "DistMax", 1.5*FMR)