Rate‐dependent fracture simulation of viscoelastic material by the phase‐field method

The phase‐field method has been studied for fracture analysis during the last decade. Based on the classical GRIFFITH‐type brittle fracture, good agreements are obtained by comparing to other numerical methods as well as experiments. Furthermore, phase‐field modelling does not depend on any explicit criterion and gets rid of discontinuous displacement tracing. Regarding most elastomers exhibiting both elastic and viscous behaviour simultaneously, it yields a rate‐dependent mechanical response. Investigated by experiments, the fracture process is shown to be rate‐dependent as well. In this work, a viscoelastic rheological model based on REESE & GOVINDJEE [1] is coupled to the phase‐field approach to investigate rate‐dependent fracture evolution within elastomers. The fracture mechanism is based on the volumetric‐isochoric split and the driving force for fracture is defined to be the elastic strain energy potential by both the equilibrium and the non‐equilibrium branches. Representative numerical models are simulated and related findings and potential subsequent perspectives are summarized to close the paper.


Introduction
Recently, phase-field modelling is considered to simulate the rate-dependent fracture evolution. A pseudo fracture resistance introduced to the strong form of the phase-field evolution, a rate-dependent fracture evolution is obtained. Meanwhile, the numerical stabilisation is also largely enhanced. Representative applications are referred to the work of MIEHE et al. [2] in quasi-static situations and the work of KUHN & MÜLLER [3] in transient situations. Moreover, motivated by experimental findings from BISCHOFF & PERRY [4] and FINEBERG et al. [5], YIN et al. [6] formulate a rate-dependent fracture toughness for the fracture evolution, obtaining good agreements with experimental data. Nevertheless, few formulations of the phasefield method to viscoelastic materials are studied. SHEN et al. [7] introduce the viscous energy dissipation as an additional driving force to evolve fracture. This method is formulated at small strain and insufficient to capture the behavior of elastomer material. SCHÄNZEL [8] formulates inelasticity in the logarithmic strain space at finite deformation using HENCKY strain, which may not define the elastic strain energy density function at the non-equilibrium branch randomly.
In order to investigate rate-dependent fracture within more general viscoelastic materials, in this work, a standard phasefield method is coupled to a viscoelastic rheological model based on REESE & GOVINDJEE [1] without considering volumetric viscosity. Due to the independent constitutive framework of the equilibrium and the non-equilibrium branches, the isochoric energy density functions for both branches can be defined flexibly, e.g. identical or not. The driving force for fracture evolution is defined by SCHÄNZEL [8], which is composed by the volumetric potential and the isochoric potential from both equilibrium and non-equilibrium branches in case of the viscoelastic solid is undergoing tension. The rheology model of the standard MAXWELL element is shown in Fig. 1. The total HELMHOLTZ free energy density function is the sum of all the elastic strain energy contributions

Constitutive model for viscoelasticity
including the volumetric energy U (J) and the isochoric free energy in the equilibrium and the non-equilibrium branches Φ eq b andΦ ne b e , respectively. The deformation gradient at the non-equilibrium branch isF e , which holdsF =F eFi . The isochoric part of the deformation gradient is defined asF = J − 1 3 F. A commonly used function for the volumetric energy is given as where κ is a model parameter known as the bulk modulus. The isochoric free energy functions for both the equilibrium and the non-equilibrium branches are based on the NEO-HOOKEAN material model, which are defined isothermally as The second axiom of thermodynamics describes the dissipation inequality of a closed system as where S represents the second PIOLA-KIRCHHOFF stress and d is the symmetric EULERIAN rate of the deformation tensor. By further extension and simplification, we yield whereτ ne iso is the KIRCHHOFF stress in the non-equilibrium branch andd i is rate of the inelastic deformation tensorF i . Based on REESE & GOVINDJEE [1], without considering the volumetric viscous contribution, the symmetric EULERIAN rate of the inelastic deformation tensord i at the non-equilibrium branch is defined as where η 0 is a model parameter describing the deviatoric viscous response. By further derivation and simplification, is obtained. Considering the alternative expressions ofb e as well asb tr ē b e = 3 a=1λ e 2 a n a ⊗ n a andb tr e = 3 a=1λ e,tr 2 a n a ⊗ n a , and take the logarithmic operation of the eigenvalues ε e a = lnλ e a and ε e,tr a = lnλ e,tr a , ε e a ≈ −η 0 devτ ne a ∆t + ε e,tr a are obtained. By iteratively solving the residual equation of a converged result of ε e a is obtained once r a ≤ T OL is reached. The detailed algorithmic setup of the local iteration is shown in [9]. Eventually, the left CAUCHY-GREEN tensor at the non-equilibrium branchb e is evaluated robustly.

Phase-field evolution
Within a continuous solid domain, a sharp crack is numerically represented by a diffusive phase-field distribution. An order parameter, described by the scalar field d(x, t), is employed to identify the material state. An unbroken solid is characterised as d = 0, while d = 1 represents the fully cracked state. The length scale parameter l is considered to govern the width of the transition zone between fractured and unfractured state. Regarding multi-dimensional problems, the crack surface energy density function reads Within the FEM application, the dependence of the diffusivity of the phase-field distribution and the length scale has been investigated by HOFACKER [10]. As a result, a minimum element size h e ≤ l/2 in the crack region regarding the finite element discretization is necessary in order to achieve realistic results. A constitutive evolution for the phase-field modelling is developed based on the energetic equilibrium. Considering a quasi-static situation within an isotropic solid domain Ω, the total internal energy density consists of the effective strain energy and the fracture energy where, G c is an intrinsic model parameter known as the critical energy release rate. In this contribution, G c is assumed to be constant and independent of other variables. The function g (d) is the degradation function to degrade the fracture driving part Φ + 0 . A variational principle is applied to the global energetic function in the reference configuration. Based on the balance of energy, an equilibrium yields The input power is P int = DW Dt . The kinetic energy of the continuum solid is defined as The rate of external work P ext is defined as The model parameter χ is defined to evaluate the viscous resistance. By simplification, the strong forms for both the mechanical response and the phase-field evolution are derived as within the volume Ω : ρ 0ü − Div g (d) P + + P − − B = 0, at the surface ∂Ω : and within the volume Ω : respectively. In order to avoid healing of the phase-field during fracture evolution, DIRICHLET boundary conditions are imposed on the phase-field degree of freedom as soon as the material is assumed to be cracked.

Numerical example
A two-dimensional boundary value problem is proposed to validate the fracture simulation of the viscoelastic solid. The setup of the geometry and the boundary conditions is clearly depicted in Fig. 2. The top and bottom edges are fully fixed along all directions and the top edge is subjected to a vertical displacement, which leads to the tensile fracture of the specimen. Two horizontal cracks are at the middle of both the left and right edges with the initial length of 20 mm. Due to the symmetrical properties, this model can be effectively simplified to a quarter of the original specimen. The material parameters are given by κ = 2.96 MPa, µ eq = 0.163 MPa, µ ne = 0.14 MPa and η 0 = 10 (s · MPa) −1 , respectively. Different loading rates are applied to investigate the rate dependent fracture evolution. The material parameters for the phase-field evolution are G c = 4.36 N/mm, l = 0.5 mm and χ = 1e −16 . The inertia effects are not taken into account in this example (ρ = 0 g/mm 3 ). The fracture evolution process is shown in Fig. 2 and the reaction forces regarding different loading rates are compared in Fig.  3. It is apparently seen that the peak loads increase with the loading rates.