Crack path comparisons of a mixed phase‐field fracture model and experiments in punctured EPDM strips

Working on quasi‐static phase‐field fracture modeling in nearly incompressible solids for crack propagation is a challenging task. To avoid arising locking effects therein, a mixed form for the solid displacement equation is developed, resulting in two unknowns: a displacement field and a hydro‐static pressure variable. In order to fulfil an inf‐sup condition, stable Taylor‐Hood elements are employed for the displacement‐pressure system. The irreversibility condition of the crack evolution is handled by help of a primal‐dual active set method. To get both a sharper crack and reasonable computational costs, adaptive meshes are used based on a predictor‐corrector scheme. The crack paths from the numerical simulations are compared on the experimentally observed crack paths in carbon black filled ethylene propylene diene monomer (EPDM) rubber strips. The punctured EPDM strips with a hole and a given notch at different heights are stretched till total failure.


Introduction
A well-established variational approach for Griffith's [1] quasi-static brittle fracture was introduced by Francfort and Marigo [2] in 1998. Since then, the method was applied in numerous different studies in calculus of variations, numerical analysis, and engineering. Phase-field modeling as introduced by Miehe et al. [3] is based on a smoothed indicator function called phasefield variable describing the crack path. If the considered solid is nearly incompressible, e.g. EPDM rubber, the quasi-static phase-field fracture model in its primal form fails due to locking effects arising in the displacement field. To allow simulating crack growth in rubber-like materials, the phase-field fracture model is extended [4,5]. The approach builds on a mixed form of the solid-displacement equation resulting in a displacement field and a pressure variable. Via Taylor-Hood elements we achieve a stable discretization of the displacement-pressure system. In this work, the developed mixed problem formulation is applied to compare the crack paths to a benchmark in EPDM rubber. The benchmark consists of punctured strips with a given notch on the left boundary which are stretched until total failure. Considering phase-field fracture modeling, not just the incompressibility, but further holes or inclusions in a material are challenging and can affect the crack path, depending on the model formulation [6], the discretization and the length scale. The main contributions of this work are: • Avoiding locking effects considering fractures in EPDM rubber by help of a mixed phase-field fracture model based on the energy functional from Wu [7] and restraining high pressure values in the inner of the crack; the energy density is split along to Amor et al. [8] (Section 2) • Solving the coupled problem with the help of a monolithic approach, a stable discretization and handling the irreversibility constraint (Section 3) • Comparing the crack paths of the numerical simulation and the experimental results in punctured EPDM strips tested with three different heights of the notch (Section 4) 2 A quasi-static phase-field model for nearly incompressible solids To allow for a weak problem formulation, we consider a subdivision 0 = t 0 < t 1 < . . . < t N = T of the interval I. In each time step, we define approximations (u n , ϕ n ) ≈ (u(t n ), ϕ(t n )) and hence the irreversibility condition is approximated by ϕ n ≤ ϕ n−1 for all n = 1, . . . , N . If it is clear from the context, we omit the index and set u := u n and ϕ := ϕ n . The phase-field space is W := H 1 (Ω) with a convex subset K := {ψ ∈ W | ψ ≤ ϕ n−1 ≤ 1}. Further, we define the function spaces V := (H 1 0 (Ω)) 2 := {w ∈ (H 1 (Ω)) 2 | w = 0 a.e. on Γ D } and U := L 2 (Ω). In the following, the critical energy release rate is denoted by G c . To guarantee well-conditioning of the system of equations, a degradation function is defined as g(ϕ) := (1 − κ)ϕ 2 + κ, with a small regularization parameter κ > 0. The stress tensor σ(u) is given by σ(u) := 2µE lin (u) + λtr(E lin (u))I with the Lamé coefficients µ, λ > 0. The linearized strain tensor therein is defined as E lin (u) := 1 2 (∇u + ∇u T ). By I, the two-dimensional identity matrix is denoted. We refer to Wu [7] for a unified phase-field fracture model with the energy functional where > 0 describes the bandwidth of the transition zone between broken and unbroken material. Based on the elastic energy functional defined in Equation (1), a mixed phase-field fracture model is presented in the following. Let a hydro-static pressure p ∈ U be defined as such that the pure elasticity equation in weak form reads as: Assume ϕ ∈ K to be given. Find u ∈ V and p ∈ U such that Remark 2.1 To avoid non-physical pressure values in the inner fracture zone, the degradation term g(ϕ) is neglected in the second term of Equation (3). So, if we are in the broken zone ϕ = 0, we are not fully divergence-free.
Along to Amor et al. [8], we consider the volumetric and deviatoric contributions of the elastic energy density separated into σ + and σ − to prevent interpenetration of the crack faces under compression [9]. For this reason, the positive part of the pressure p + ∈ L 2 (Ω) as well as the positive part of E + lin (u) has to be defined as p + := max{p, 0} and E + lin (u) := max{E lin (u), 0} as the maximum function, such that the stress tensor σ(u, p) is split into: Next, we can formulate the mixed phase-field problem in incremental form similar to [4]: Problem 2.2 (Mixed Phase-field Formulation) Let the initial data ϕ n−1 ∈ K be given. Find u := u n ∈ {u D + V}, p := p n ∈ U and ϕ := ϕ n ∈ K for the loading steps n = 1, 2, . . . , N such that

Remark 2.3
In the elasticity part, time-lagging is used in the phase-field variable ϕ to obtain a convex functional, see [10] for further details.
The numerical treatment of the phase-field system in Problem 2.2 in a monolithic fashion including the discretization as well as the adaptive refinement strategy are discussed in the following.

Solving and implementation
The crack irreversibility is respected with a primal-dual active-set method, see [11] for further details. For the spatial discretization, we employ a Galerkin finite element method in each incremental step, where the domain Ω is partitioned into quadrilaterals. To fulfill a discrete inf-sup condition, stable Taylor-Hood elements with biquadratic shape functions (Q 2 ) for the displacement field u and bilinear shape functions (Q 1 ) for the pressure variable p and the phase-field variable ϕ are used as in [4]. The open-source code to derive the proposed problem formulation is available at https://github.com/tjhei/cracks which is embedded in deal.II [12]. Further details on the code are given in [13]. For adaptive refined meshes, a predictor-corrector scheme is used along to [11,Chapter 4]. This refinement scheme allows to refine the mesh locally depending on a propagating fracture with a choosen threshold for the phase-field variable.

Comparing numerical and experimental results in punctured EPDM strips 4.1 EPDM benchmark
The experiments are conducted with strips of Sulphur crosslinked EPDM (Keltan 2450) filled with 60 phr carbon black N550. The strips have a given notch on the left side and a hole with a diameter of 8mm. Crack propagation is observed via stretching the punctured strips until total failure with a speed of 200mm/min. To allow simulating the crack propagation in punctured strips in a two-dimensional setup by help of the phase-field fracture model from Section 2, we reduce the geometry to the area of interest between the bulges on the bottom and top part where the specimens are fixed. The geometry for the numerical simulation is depicted in Figure 1. We choose homogeneous Dirichlet boundary conditions u y = 0 on the top boundary and fix the strips in horizontal direction on the top and bottom boundary via u x = 0. Further, the following boundary condition characterizes the loading forces on the bottom boundary Γ force in y-direction: where t denotes the total time. The time interval I is divided into steps of size δt. The given notch is discretized with an initial condition ϕ = 0. Further material and numerical parameters are listed in Table 1.    The material parameters λ and µ as well as G c are identified via evaluating additional experiments. Hereby, G c is identified via a single edge notched tension test in pure shear mode using "Digital Image Correlation" via GOM Correlate. In case of the first loading of this specific EPDM material non-linearity up to 100% strain plays a minor role. Although geometrical non-linearity is relevant in general in terms of elastomers, it is neglected in the model for now.

Numerical results
Based on the setup of the EPDM benchmark from the previous section, the numerical results of three test cases are presented. We consider a given notch at three different heights from the bottom line of the EPDM strips: 6, 12, and 18mm. With a height of 18mm the notch is exactly on the height of the center of the inclusion. In all tests, using the predictor-corrector scheme from [11], three adaptive refinement cycles per incremental step are allowed with a phase-field threshold of 0.7. In Figure 2 with a given notch at 6mm it is observed, that the inclusion in the material has an impact on the direction of the crack path in the experiment as well as in the numerical simulation.
With a given notch at 12mm as depicted in Figure 3, the crack propagates to the inclusion and leaves the hole approximately straight. In the Figures 3 and 4, the experimental results show different behaviour in the crack path from the hole to the right boundary. One reason can be inhomogeneities in the EPDM rubber material. www.gamm-proceedings.com

Conclusion
In this work, a phase-field fracture model in mixed form is proposed to simulate crack propagation in EPDM rubber considering the crack path behavior. The model is based on a mixed form of the solid equation derived from Wu's energy functional and the approach of Amor et al. for a volumetric-deviatoric stress splitting matching the volume conservation of incompressible materials. The obtained crack path in the numerical results compared to the punctured EPDM strips after stretching them with a certain load till total failure, coincide aside from experimental scattering. In future work, further quantities of interest such as load-displacement curves, crack length and crack and bulk energy will be evaluated to validate the mixed phase-field fracture model in detail. Further, to improve the efficiency of the numerical solving, we are on the way to develop a preconditioner.