Ductile fracture modeling by phase field, Hencky strain elasticity and finite J2 plasticity using nonlocal operator method

TL;DR

This paper introduces a novel phase field model for ductile fracture that integrates Hencky strain, finite J2 plasticity, and the nonlocal operator method, providing a comprehensive approach to simulate ductile failure.

Contribution

It develops a variational phase field model incorporating advanced plasticity and strain measures with a spectral decomposition algorithm, enhancing simulation accuracy for ductile fracture.

Findings

Validated with numerical examples including notched plates and necking bars.

Demonstrated the effectiveness of the model in capturing ductile fracture behavior.

Provided a consistent computational framework with derived stiffness matrices.

Abstract

A phase field model for ductile fracture considering Hencky strain and finite J2 plasticity is presented using the nonlocal operator method. A variational derivation of J2 plasticity at finite strain with a phase field model is performed. The method includes a logarithmic strain tensor and an exponential mapping in the plasticity evolution. A spectral decomposition based algorithm for computing the first and second order derivatives of the composite matrix function is implemented. A consistent tangential stiffness matrix is derived and used in Newton-Raphson iterations. Several numerical examples are performed to validate the method, including notched single-edged plates with brittle fracture or ductile fracture and necking of a bar with/without phase field model.

Tables

Table 1. Table 1: Necking of a rod

Material parameters
Bulk modulus	$\kappa =164.206$ kN/mm2
Shear modulus	$\mu =80.1938$ kN/mm2
Initial yield stress	$y_0=0.45$ kN/mm2
Infinite yield stress	$y_{\infty }=0.715$ kN/mm2
Hardening modulus	$h=0.12924$ kN/mm2
Saturation parameter	$\Delta =16.93$

Equations

$$\underset{u,K}{\min }\left\{\alpha {\int }_{\Omega }|u-\bar{u}{|}^2+{\int }_{\Omega \backslash K}|\nabla u{|}^2+\beta {\mathcal{H}}^{n-1}(K)\right\},$$

$$A{T}_{ϵ}(u,c)=$$

$${\mathbf{K}}_i:={\int }_{{\mathcal{S}}_i}w({\mathbf{r}}_{ij}){\mathbf{r}}_{ij}\otimes {\mathbf{r}}_{ij}d{V}_j,$$

$${\mathcal{S}}_i^{\prime }=\{{\mathbf{X}}_j|{\mathbf{X}}_i\in {\mathcal{S}}_j\}.$$

$$\tilde{\nabla }\ast {\mathbf{u}}_i:={\int }_{{\mathcal{S}}_i}\omega (r_{ij})({\mathbf{K}}_i^{-1}{\mathbf{r}}_{ij})\ast {\mathbf{u}}_{ij}d{V}_j,$$

$$\Omega =\sum _{i=1}^N{V}_i,$$

$${\mathcal{N}}_i=\{i,j_1,..,j_k,..,j_{n_i}\},$$

$$\tilde{\nabla }{u}_i={\int }_{{\mathcal{S}}_i}\omega (r_{ij}){\mathbf{K}}_i^{-1}{\mathbf{r}}_{ij}(u_j-u_i)d{V}_j=\sum _{j\in {\mathcal{N}}_i}{V}_j\omega (r_{ij}){\mathbf{K}}_i^{-1}{\mathbf{r}}_{ij}(u_j-u_i).$$

$$V_j\omega (r_{ij}){\mathbf{K}}_i^{-1}{\mathbf{r}}_{ij}=\{\begin{array}{l}(r_j^x,r_j^y)\text{mbox}in2D\\ (r_j^x,r_j^y,r_j^z)\text{mbox}in3D.\end{array}$$

$$\tilde{\nabla }{u}_i=\left(\begin{array}{c}\frac{\partial u_i}{\partial x}\\ \frac{\partial u_i}{\partial y}\\ \frac{\partial u_i}{\partial z}\end{array}\right)=\underset{{\mathbf{B}}_i}{\underbrace{\left(\begin{array}{cccc}-{\sum }_{j\in {\mathcal{N}}_i}{r}_j^x& r_{j_1}^x& \dots & r_{j_n}^x\\ -{\sum }_{j\in {\mathcal{N}}_i}{r}_j^y& r_{j_1}^y& \dots & r_{j_n}^y\\ -{\sum }_{j\in {\mathcal{N}}_i}{r}_j^z& r_{j_1}^z& \dots & r_{j_n}^z\end{array}\right)}}\cdot \underset{{\mathbf{U}}_i}{\underbrace{\left(\begin{array}{c}{u}_i\\ u_{j_1}\\ ⋮\\ u_{j_n}\end{array}\right)}}.$$

$${\mathcal{F}}^{hg}$$

$$=\frac{1}{2}\frac{\alpha }{\text{mbox}tr({\mathbf{K}}_i)}\sum _{j\in {\mathcal{N}}_i}\omega (r_{ij})(\tilde{\nabla }{u}_i{\mathbf{r}}_{ij}-u_{ij})\cdot (\tilde{\nabla }{u}_i{\mathbf{r}}_{ij}-u_{ij})V_j,$$

$${\mathcal{F}}^{hg}$$

$$\displaystyle\mathbb{K}_{i}^{hg}:=\frac{\alpha}{\mbox{tr}(\bm{K}_{i})}\Big{(}\begin{pmatrix}\sum_{j\in\mathcal{N}_{i}}\omega(r_{ij})V_{j}&-\omega(r_{ij_{1}})V_{j_{1}}&\cdots&-\omega(r_{ij_{n}})V_{j_{n}}\\ -\omega(r_{ij_{1}})V_{j_{1}}&\omega(r_{ij_{1}})V_{j_{1}}&\bm{0}&0\\ \vdots&\bm{0}&\ddots&\bm{0}\\ -\omega(r_{ij_{n}})V_{j_{n}}&0&\bm{0}&\omega(r_{ij_{n}})V_{j_{n}}\end{pmatrix}-\bm{B}_{i}^{T}\bm{K}_{i}\bm{B}_{i}\Big{)}.$$

$$\left(\begin{array}{c}\frac{\partial u_i}{\partial x}\\ \frac{\partial u_i}{\partial y}\\ \frac{\partial u_i}{\partial z}\\ \frac{\partial v_i}{\partial x}\\ \frac{\partial v_i}{\partial y}\\ \frac{\partial v_i}{\partial z}\\ \frac{\partial w_i}{\partial x}\\ \frac{\partial w_i}{\partial y}\\ \frac{\partial w_i}{\partial z}\end{array}\right)=\underset{{\mathbf{B}}_i}{\underbrace{\left(\begin{array}{c}\text{-}{\sum }_{j\in {\mathcal{N}}_i}{r}_j^x,0,0,r_{j_1}^x,0,0,\cdots ,r_{j_n}^x,0,0\\ \text{-}{\sum }_{j\in {\mathcal{N}}_i}{r}_j^y,0,0,r_{j_1}^y,0,0,\cdots ,r_{j_n}^y,0,0\\ \text{-}{\sum }_{j\in {\mathcal{N}}_i}{r}_j^z,0,0,r_{j_1}^z,0,0,\cdots ,r_{j_n}^z,0,0\\ 0,\text{-}{\sum }_{j\in {\mathcal{N}}_i}{r}_j^x,0,0,r_{j_1}^x,0,\cdots ,0,r_{j_n}^x,0\\ 0,\text{-}{\sum }_{j\in {\mathcal{N}}_i}{r}_j^y,0,0,r_{j_1}^y,0,\cdots ,0,r_{j_n}^y,0\\ 0,\text{-}{\sum }_{j\in {\mathcal{N}}_i}{r}_j^z,0,0,r_{j_1}^z,0,\cdots ,0,r_{j_n}^z,0\\ 0,0,\text{-}{\sum }_{j\in {\mathcal{N}}_i}{r}_j^x,0,0,r_{j_1}^x,\cdots ,0,0,r_{j_n}^x\\ 0,0,\text{-}{\sum }_{j\in {\mathcal{N}}_i}{r}_j^y,0,0,r_{j_1}^y,\cdots ,0,0,r_{j_n}^y\\ 0,0,\text{-}{\sum }_{j\in {\mathcal{N}}_i}{r}_j^z,0,0,r_{j_1}^z,\cdots ,0,0,r_{j_n}^z\end{array}\right)}}\underset{{\mathbf{U}}_i}{\underbrace{\left(\begin{array}{c}{u}_i\\ v_i\\ w_i\\ u_{j_1}\\ v_{j_1}\\ w_{j_1}\\ ⋮\\ v_{j_n}\\ w_{j_n}\end{array}\right)}}.$$

$${\mathcal{F}}^{hg}$$

$$=\frac{1}{2}\frac{\alpha }{\text{mbox}tr({\mathbf{K}}_i)}\sum _{j\in {\mathcal{N}}_i}\omega (r_{ij})(\tilde{\nabla }{\mathbf{u}}_i{\mathbf{r}}_{ij}-{\mathbf{u}}_{ij})\cdot (\tilde{\nabla }{\mathbf{u}}_i{\mathbf{r}}_{ij}-{\mathbf{u}}_{ij})V_j,$$

$${\mathcal{F}}^{hg}$$

$$\displaystyle\mathbb{K}_{i}^{hg}:=\frac{\alpha}{\mbox{tr}(\bm{K}_{i})}\Big{(}\begin{pmatrix}\sum_{j\in\mathcal{N}_{i}}\omega(r_{ij})V_{j}&\text{-}\omega(r_{ij_{1}})V_{j_{1}}&\cdots&\text{-}\omega(r_{ij_{n}})V_{j_{n}}\\ \text{-}\omega(r_{ij_{1}})V_{j_{1}}&\omega(r_{ij_{1}})V_{j_{1}}&\bm{0}&0\\ \vdots&\bm{0}&\ddots&\bm{0}\\ \text{-}\omega(r_{ij_{n}})V_{j_{n}}&0&\bm{0}&\omega(r_{ij_{n}})V_{j_{n}}\end{pmatrix}\otimes\bm{I}-\bm{B}_{i}^{T}(\bm{K}_{i}\otimes\bm{I})\bm{B}_{i}\Big{)}.$$

$$e^x=\sum _{n=0}^{\infty }\frac{x^n}{n!},e^{\mathbf{X}}:=\sum _{n=0}^{\infty }\frac{{\mathbf{X}}^n}{n!}.$$

$$\frac{dy}{dt}=ay\to y(t)=e^{at}y(0)$$

$$\frac{d\mathbf{y}}{dt}=\mathbf{A}\mathbf{y}\to \mathbf{y}(t)=e^{\mathbf{A}t}\mathbf{y}(0),$$

$$\mathbf{A}=\left(\begin{array}{ccc}{a}_{11}& a_{12}& a_{13}\\ a_{12}& a_{22}& a_{23}\\ a_{13}& a_{23}& a_{33}\end{array}\right),$$

$$\mathbf{A}=\sum _{i=1}^3{\lambda }_i{\mathbf{n}}_i\otimes {\mathbf{n}}_i,$$

$$f(\mathbf{A})=\sum _{k=1}^{3}f({\lambda }_k){\mathbf{n}}_k\otimes {\mathbf{n}}_k=:\left(\begin{array}{ccc}{f}_{11}& f_{12}& f_{13}\\ f_{12}& f_{22}& f_{23}\\ f_{13}& f_{23}& f_{33}\end{array}\right),$$

$$\displaystyle\frac{\partial f(\bm{A})}{\partial a_{ij}}=\sum_{k=1}^{3}\Big{(}\frac{\partial f(\lambda_{k})}{\partial a_{ij}}\bm{n}_{k}\otimes\bm{n}_{k}+f(\lambda_{k})\big{(}\frac{\partial\bm{n}_{k}}{\partial a_{ij}}\otimes\bm{n}_{k}+\bm{n}_{k}\otimes\frac{\partial\bm{n}_{k}}{\partial a_{ij}}\big{)}\Big{)}.$$

$$\left(\begin{array}{ccc}\frac{\partial f_{11}}{\partial a_{11}}& \frac{\partial f_{11}}{\partial a_{12}}& \frac{\partial f_{11}}{\partial a_{22}}\\ \frac{\partial f_{12}}{\partial a_{11}}& \frac{\partial f_{12}}{\partial a_{12}}& \frac{\partial f_{12}}{\partial a_{22}}\\ \frac{\partial f_{22}}{\partial a_{11}}& \frac{\partial f_{22}}{\partial a_{12}}& \frac{\partial f_{22}}{\partial a_{22}}\end{array}\right)=\left(\begin{array}{ccc}{f}^{\prime }(\lambda )& 0& 0\\ 0& f^{\prime }(\lambda )& 0\\ 0& 0& f^{\prime }(\lambda )\end{array}\right).$$

$$\partial {\lambda }_i={\mathbf{n}}_i^T\partial (\mathbf{A}){\mathbf{n}}_i,\partial {\mathbf{n}}_i={\left({\lambda }_i\mathbf{I}-\mathbf{A}\right)}^{+}\partial (\mathbf{A}){\mathbf{n}}_i,$$

$$\frac{\partial {\lambda }_i}{\partial \mathbf{A}}={\mathbf{n}}_i\otimes {\mathbf{n}}_i,\frac{\partial {\mathbf{n}}_i}{\partial \mathbf{A}}=\sum _{j\ne i}^3\frac{1}{{\lambda }_i-{\lambda }_j}{\mathbf{n}}_j\otimes {\mathbf{n}}_j\otimes {\mathbf{n}}_i$$

$$\frac{\partial f(\mathbf{A})}{\partial b_l}=\frac{\partial f(\mathbf{A})}{\partial a_{ij}}\frac{\partial a_{ij}}{\partial b_l},$$

$$\displaystyle\frac{\partial^{2}f(\bm{A})}{\partial a_{ij}\partial a_{mn}}=\sum_{k=1}^{3}\Big{(}\frac{\partial^{2}f(\lambda_{k})}{\partial a_{ij}\partial a_{mn}}\bm{n}_{k}\otimes\bm{n}_{k}+\frac{\partial f(\lambda_{k})}{\partial a_{ij}}\frac{\partial\bm{n}_{k}}{\partial a_{mn}}\otimes\bm{n}_{k}$$

$$\displaystyle+\frac{\partial f(\lambda_{k})}{\partial a_{ij}}\bm{n}_{k}\otimes\frac{\partial\bm{n}_{k}}{\partial a_{mn}}+\frac{\partial f(\lambda_{k})}{\partial a_{mn}}\big{(}\frac{\partial\bm{n}_{k}}{\partial a_{ij}}\otimes\bm{n}_{k}+\bm{n}_{k}\otimes\frac{\partial\bm{n}_{k}}{\partial a_{ij}}\big{)}$$

Full text

1 Introduction

Safety evaluation of metallic materials or structures is important for many engineering applications. In the process of determining the critical state of the structures, the formation of plastic deformation and the development of brittle/ductile fractures are often involved. Plastic deformation is the ability of a solid material to permanently deform under the action of external forces when the critical state of elastic deformation reaches yielding. Ductile fracture with finite deformation means that the fracture develops in conjunction with the permanent deformation and the possible thermal effect. The simulation of ductile fracture encounters two difficulties: the finite elastoplastic deformation and the evolution of the fractures.

The mechanism of elastoplastic deformation is closely related to micromechanics based on the internal structure of materials. However, a phenomenological theory of plasticity based on a macroscopic description is useful for many engineering applications. Finite strain plasticity is usually based on the multiplicative decomposition of the deformation gradient proposed by Lee [32, 33]. The multiplicative decomposition consists of an elastic and a plastic part of the deformation tensor and allows a material frame invariant description of plasticity [27]. The thermodynamically consistent finite strain plasticity framework has been well presented by Havner, Hill, Rice, and others [21, 56, 6, 45, 17].

The plastic model mainly considers plastic deformation, and the evolution of sharp discontinuities such as fractures requires additional techniques such as the phase field model. Phase field methods have proven to be very successful in modeling brittle fractures and ductile fractures. Based on the Griffith energy principle, the fracture surface is represented by a phase field energy functional in terms of $\Gamma$-convergence. The evolution of the phase field reflects the automatic determination of the fracture direction and evolution. The development of the phase field energy framework took advantage of the Mumford-Shah functional [44] and the Ambrosio-Tortorelli functional [3]. Mumford-Shah functional was introduced in [44] for image segmentation, which is written as

[TABLE]

where $u\in\Omega\backslash K$, $\bar{u}$ the prescribed field value of $u$, $K\subset\Omega$ closed in $\Omega,\mathcal{H}^{n-1}(K)<\infty$, where $\mathcal{H}^{n-1}$ is the ($n-1$)-dimensional Hausdorff measure (e.g. surface area in 3D when $n=3$), and $\alpha,\beta>0$ are fixed parameters. The functional consists of the sharp interface in domain $K$ and low gradient terms in $\Omega\backslash K$. In Ref [3], Ambrosio and Tortorelli proposed an elliptic functional (i.e. the Ambrosio-Tortorelli functional) as follows:

[TABLE]

where $\epsilon$ is a positive fixed parameter. The authors proved that $AT_{\epsilon}(u,c)$ converges in the sense of $\Gamma$-convergence with the Mumford-Shah functional when $\epsilon\to 0$. In the spirit of the Ambrosio-Tortorelli functional, Francfort and Marigo [23] replaced the gradient term $|\nabla u|^{2}$ in Eq.2 by the elastic strain energy density in solid mechanics and formulated the brittle fracture as an energy minimization problem. The numerical experiments were performed in [10]. The additional field $c$ greatly simplifies the description of the fracture. In Ref [43], Miehe etal developed the thermodynamically consistent phase field model based on the spectral decomposition of the strain tensor in isotropic elastic material, showing great stability and capability in complicated crack simulations. Phase field models have been used to solve brittle fractures [29, 2], hydraulic fractures [58, 41, 61, 26], ductile fractures with plasticity [8, 7, 19, 25, 1], geomaterials [16, 22, 12, 60], plate and shell [5, 48], multiphysics problems [42, 39, 18] and many others. For a more comprehensive review of the phase field model, the reader is referred to [59, 13, 62].

Numerical methods for the phase field model are usually based on the finite element method [43] or isogeometric analysis [9, 15]. Nonlocal Operator Method (NOM) is a numerical method based on nonlocality to solve partial differential equations [49, 53, 54, 52]. NOM as a generalization of the dual-horizon Peridynamics (PD) [50], extends the first-order nonlocal gradient to the higher-order nonlocal derivatives without recourse to shape functions. By adopting a variational framework and weighted residual methods, NOM is applicable to many physical problems based on the weak form [55]. The implementation of NOM allows an inhomogeneous discretization and has no restriction on the support size and support shape. Meanwhile, it has the capability of deriving the nonlocal models based on the traditional local forms or the weak forms [51]. It has been applied to solve various physical problems such as linear/nonlinear gradient elasticity [54], waveguide problem in electromagnetic field [49], Cahn-Hilliard equation [55], von-Karmman equations for thin plates [52]. In this work, we apply the NOM to model ductile fracture due to plasticity with consideration of both material and geometric nonlinearity.

The remainder of the paper is organized as follows. In Section 2, the first- order NOM is briefly discussed and the matrix form, the tangent stiffness matrix of the hourglass energy, are explicitly formulated with a discrete setting of the support domain. In Section 3, the spectral decomposition of the matrix function is presented and its first-order and second- order derivatives are discussed. In Section 4, the governing equation of Hencky strain elasticity is derived and the residual vector and tangent stiffness matrix are derived. In Section 5, the theoretical framework of multiplicative J2 elastoplasticity including flowing rule and yield function is presented. In Section 6, the phase field J2 plasticity model is discussed. Section 7 presents several numerical examples, including a notched-single-edge plate with brittle fracture or ductile fracture and the necking of a bar with/without ductile fracture, to validate the NOM scheme. Some conclusions are drawn in Section 8.

2 Brief review of nonlocal operator method

2.1 Support, dual-support and nonlocal operators

Consider a domain as shown in Fig.1, let $\bm{X}$ be spatial coordinates in the initial configuration $\mathbf{\Omega}$; $\bm{r}_{ij}:=\bm{X}_{j}-\bm{X}_{i}$ is a vector (or a spatial vector, or simply a vector) starts from $\bm{X}_{i}$ to $\bm{X}_{j}$; $\bm{v}_{i}:=\bm{v}(\bm{X}_{i},t)$ and $\bm{v}_{j}:=\bm{v}(\bm{X}_{j},t)$ are the field values for $\bm{X}_{i}$ and $\bm{X}_{j}$, respectively; $\bm{v}_{ij}:=\bm{v}_{j}-\bm{v}_{i}$ is the relative field vector for spatial vector $\bm{r}_{ij}$.

Support $\mathcal{S}_{i}$ is a finite-size neighborhood of point $\bm{X}_{i}$. Support $\mathcal{S}_{i}$ specifies the range of nonlocal interaction that happened with respect to point $\bm{X}_{i}$. The main function of a support is to define different nonlocal operators. In mathematics, a specific quantitative measure of the support can be described by a moment (or a shape tensor),

[TABLE]

where $w(\bm{r}_{ij})$ is the weight function.

Dual-support is defined as a union of the points whose supports include $\bm{X}_{i}$, denoted by

[TABLE]

The point $\bm{X}_{j}$ forms dual-vector $\bm{r}_{ij}^{\prime}(=\bm{X}_{i}-\bm{X}_{j}=-\bm{r}_{ij})$ in $\mathcal{S}_{i}^{\prime}$. On the other hand, $\bm{r}_{ij}^{\prime}$ is the spatial vector formed in $\mathcal{S}_{j}$. One example to illustrate the range of support and dual-support is shown in Fig.1. Seven representative points and their support areas are drawn in the figure. After defining the support area, we can see that the support of $\bm{X}$ contains $\{\bm{X}_{1},\bm{X}_{2},\bm{X}_{4},\bm{X}_{6}\}$, while the dual-support of $\bm{X}$ consists of $\{\bm{X}_{1},\bm{X}_{2},\bm{X}_{3},\bm{X}_{4}\}$.

2.2 Nonlocal operators in support

The common operators in calculus include the gradient of the scalar and vector field, the curl and divergence of the vector field. The definitions of some nonlocal operators can be found in reference [20]. These operators have the corresponding nonlocal forms based on the Taylor series expansion. We use ${\tilde{\nabla}}$ to denote the nonlocal operator, while the local operators follow the conventional notations.

The first-order nonlocal operator for vector-valued field $\bm{u}$ at point $\bm{x}_{i}$ is defined as

[TABLE]

where $*\in\{\cdot,\otimes,\times\}$ denote the inner product, tensor product and cross product, respectively. The continuous form is inconvenient for implementation. After discretization of the domain by particles, the whole domain is represented by

[TABLE]

where $i$ is the global index of volume $V_{i}$, $N$ is the number of particles in $\Omega$.

Particles in $\mathcal{S}_{i}$ are represented by

[TABLE]

where $j_{1},..,j_{k},..,j_{n}$ are the global indices of neighbors of particle $i$. Herein, the number of particles in each support is fixed as $n$, although $n$ can vary for each particle.

The nonlocal gradient for scalar-valued field $u$ at point $\bm{X}_{i}$ in continuous form and discrete form is

[TABLE]

Let

[TABLE]

Then $\tilde{\nabla}u_{i}$ can be written as matrix form as

[TABLE]

The matrix form of $\tilde{\nabla}u_{i}$ is similar to the gradient of the shape function in the finite element method.

The energy functional for scalar-valued field $u$ at point $\bm{X}_{i}$

[TABLE]

where $\alpha$ is a penalty coefficient. $\mbox{tr}(\square)$ denotes the trace of the matrix.

Replacing $\nabla u_{i}$ with $\tilde{\nabla}u_{i}$ in Eq.8, $\mathcal{F}^{hg}$ in discrete form can be simplified as

[TABLE]

where

[TABLE]

Actually, $\mathbb{K}_{i}^{hg}$ is the hourglass matrix that can suppress the hourglass mode when adding to the tangent stiffness matrix of the physical model. More details can be found in [53].

Eq.13 is for a single field $u$. For vector-valued field $\bm{u}=(u,v,w)$ in 3-dimensional space, the nonlocal derivatives in matrix form are

[TABLE]

The corresponding energy functional at point $\bm{X}_{i}$

[TABLE]

Replacing $\nabla\bm{u}_{i}$ with $\tilde{\nabla}\bm{u}_{i}$ in Eq.5, $\mathcal{F}^{hg}$ in discrete form can be simplified as

[TABLE]

where $\bm{U}_{i}^{T}=(u_{i},v_{i},w_{i},u_{j_{1}},v_{j_{1}},w_{j_{1}},\cdots,u_{j_{n}},v_{j_{n}},w_{j_{n}})$, and $\mathbb{K}_{i}^{hg}$ is the hourglass tangent stiffness matrix for vector field $\bm{u}=(u,v,w)$:

[TABLE]

Herein, only the first order nonlocal gradient is provided. For higher order formulation of NOM, the interested reader is referred to [52].

3 Matrix function and spectral decomposition

A difficulty in single crystal plasticity with large strains is the calculation of exponential matrix functions and their derivatives for the incompressibility of the plasticity evolution. The use of exponential matrix function allows the incompressibility of plasticity evolution if the exponential matrix function is calculated accurately. There are several methods of calculating matrix functions, for example, direct summation of (truncated) series expansion of exponential matrix [46], spectral decomposition [14], generating function of tensor functions [34, 31] and so on. Lu [34] used stem function to compute the arbitrary tensor functions and their first derivatives with closed-form, singularity-free expressions. The spectral decomposition [46, 14] can provide a closed-form representation of both the exponential and its derivative. In the case of repeated eigenvalues, special care should be taken to avoid numerical singularities. The spectral decomposition has been used in several implementations of plasticity at finite-strain plasticity [28, 46]. Korelc and Stupkiewicz [31] have proposed a closed-form matrix function using automatic differentiation of an appropriate scalar generating function. High numerical accuracy is achieved even for the repeated eigenvalues. In the present work, we use spectral decomposition to compute the composite matrix function and its first-order/ second- order derivatives. An explicit simple expression is provided for the case of triply repeated eigenvalues.

3.1 Matrix function

The matrix function can be obtained by generalizing the scalar-valued function to the matrix-valued function. For example, using series expansion, the exponential matrix function can be defined as

[TABLE]

Matrix function has profound influence on vector-valued differential equations. The generalization of scalar-valued ordinary differential equation to vector-valued ordinary differential equations requires the matrix functions, for example,

[TABLE]

where $a,\bm{A}$ are coefficients, $\bm{y}\in\mathbb{C}^{n},\bm{A}\in\mathbb{C}^{n\times n}$.

3.2 Spectral decomposition for matrix function and its derivatives

In continuum mechanics, the third-order matrix is often used. Let us focus on the symmetric matrix function of third-order

[TABLE]

where $a_{ij}:=a_{ij}(b_{1},...,b_{n})$ are the functions of variables ($b_{1},b_{2},...,b_{n})$.

The eigenvalue decomposition of $\bm{A}$ is

[TABLE]

where $\lambda_{i}$ and $\bm{n}_{i}$ are the distinct eigenvalues and orthogonal unit eigenvectors, respectively. For the case of repeated eigenvalues, the distinct eigenvectors can be achieved by perturbation of the repeated eigenvalues.

Based on spectral decomposition, the matrix function of $\bm{A}$ can be written as

[TABLE]

where $f_{ij}$ are the items of the matrix function.

Using chain-rule, the derivative of $f(\bm{A})$ on ($a_{11},a_{12},a_{13},a_{22},a_{23},a_{33}$)

[TABLE]

When the eigenvalues are identical, the matrix becomes a identity matrix, the matrix derivative has a simple form

[TABLE]

3.3 Matrix derivative

Derivatives of eigenvalue and eigenvector with respect to matrix $A$ are

[TABLE]

where $\bm{A}$ is real and symmetric, $\lambda_{i}$ and $\bm{n}_{i}\left(\bm{n}_{i}\cdot\bm{n}_{j}=\delta_{ij}\right)$ are the distinct eigenvalues and eigenvectors of $\bm{A}$. The generalized inverse of ($\lambda_{i}\bm{I}-\bm{A}$) has the explicit form: $\left(\lambda_{i}\bm{I}-\bm{A}\right)^{+}=\sum_{j\neq i}^{3}\frac{1}{\lambda_{i}-\lambda_{j}}\bm{n}_{j}\otimes\bm{n}_{j}$. Then the matrix derivative can be written as well as

[TABLE]

Therefore, the derivative of matrix function $f(\bm{A})$ with respect to $b_{l}$ can be written as

[TABLE]

where $\frac{\partial f(\bm{A})}{\partial a_{ij}}$ is computed by Eq.21.

For many real applications, it is also required to calculate the second-order derivative of the matrix function with respect to $\bm{A}$, that is

[TABLE]

In the above equation, the second order derivative of eigenvalue and eigenvector with respect to the variables are required. The second order derivative can be obtained by calculating the partial derivatives of both sides of Eq.23, which is written as

[TABLE]

The above two equations contain all terms such as $\frac{\partial^{2}\lambda_{k}}{\partial a_{ij}\partial a_{mn}}$ and $\frac{\partial^{2}\bm{n}_{k}}{\partial a_{ij}\partial a_{mn}}$, which can be extracted by mathematical software such as Mathematica.

Then the second derivative of matrix function $f(\bm{A})$ with respect to $b_{l}$ can be written as

[TABLE]

4 Hyperelasticity based on Hencky strain

Hencky strain is suitable to describe the moderate large deformation of elastic solid [4]. Hencky strain is described by the logarithmic function of Cauchy-Green tensor. Using conventional small deformation material tensor, the finite strain energy density based on logarithmic strain of left Cauchy-Green tensor ($\bm{b}=\bm{F}\bm{F}^{T}$) can be written as

[TABLE]

where $\bm{\varepsilon}=\frac{1}{2}\ln\bm{b}$ is the logarithmic strain, $\bm{\sigma}$ is the stress tensor via the conventional 4th-order material tensor $\mathbb{C}$ as follows.

[TABLE]

The variation of strain energy in domain $\Omega$ is derived as

[TABLE]

where

[TABLE]

In the derivation, $\frac{\partial\phi}{\partial\bm{b}}$ is symmetric and the relation of $A:(BC)=\mbox{tr}[A^{T}BC]=B^{T}A:C=A^{T}B:C^{T}$ is used. In the derivation, the boundary term is neglected for the sake of conciseness. For any $\delta\bm{u}$, $\delta\Psi=0$ in domain $\Omega$ yields the governing equations

[TABLE]

where $\bm{g}$ is the body force density.

In order to derive the residual vector and tangent stiffness matrix of the energy functional at a point, the Voigt notation of $\bm{b}$ is used

[TABLE]

Consider energy functional $\phi:=\phi(\bar{\bm{b}}(\bm{F}))$. Consider the variation of $\phi$

[TABLE]

where

[TABLE]

The calculation of $\frac{\partial\ln\bm{b}}{\partial\bar{\bm{b}}},\frac{\partial^{2}\ln\bm{b}}{\partial\bar{\bm{b}}\partial\bar{\bm{b}}}$ should be done for each sub-term based on the spectral decomposition and be constructed based on the Voigt notation of $\bm{b}$.

So the residual vector and tangent stiffness matrix at a point can be written as

[TABLE]

where $\bm{B}$ is defined in Eq.14. In above equations, $\frac{\partial\bar{\bm{b}}}{\partial\bm{F}}$ and $\frac{\partial^{2}\bar{\bm{b}}}{\partial\bm{F}\partial\bm{F}}$ are independent of the material constitutive. When the residual vector and tangent stiffness matrix are derived and the appropriate boundary conditions are enforced, the Newton iteration scheme can be applied to find the solution.

5 Variational derivation of finite strain elastoplasticity with phase field model

5.1 Multiplicative elastoplasticity

The solid material is assumed as isotropic elasticity, whose strain energy is expressed using the invariants of the Hencky strain of the form

[TABLE]

where $\bm{b}_{e}$ is the elastic part of the left Lagrange-Green strain tensor.

In multiplicative elastoplasticity theory [57], the deformation gradient $\bm{F}$ can be decomposed into as the elastic part $\bm{F}_{e}$ and plastic part $\bm{F}_{p}$ as,

[TABLE]

Then the following standard kinematic quantities (total, elastic, and plastic Cauchy-Green tensors) are

[TABLE]

The velocity gradient, its elastic and plastic parts, are defined by

[TABLE]

The rate of deformation tensor $\bm{d}$ and spin tensor $\bm{w}$ are defined for $\bm{l}$

[TABLE]

Similarly, $\bm{l}_{p}$ and $\bm{L}_{p}$ can be decomposed into symmetric and antisymmetric parts, namely

[TABLE]

Based on the above definitions, $\dot{\overline{\bm{C}_{p}^{-1}}}$ satisfies the relationship

[TABLE]

The yield function and plastic flow rule can be defined either in the intermediate configuration or in the current configuration. For the former case, the adequate stress is the Mandel stress tensor [35, 36]

[TABLE]

where $\bm{\tau}=2\bm{b}_{e}\frac{\partial\phi_{e}}{\partial\bm{b}_{e}}$ is the Kirchhoff stress. The Mandle stress is symmetric in isotropic elasticity. The yield function of $J_{2}$ plasticity based on the Mandel stress has the form

[TABLE]

where $\bm{M}^{\prime}=\bm{M}-1/3\mbox{tr}(\bm{M})\bm{I}$, $||\bm{A}||=\sqrt{\bm{A}:\bm{A}}$ and $\sigma_{y}(\alpha)$ is the yield stress function of hardening variable $\alpha$. $\bm{\tau}_{dev}=\bm{\tau}-1/3\mbox{tr}(\bm{\tau})\bm{I}$. The associated plastic flow rule,

[TABLE]

Apparently, $\bm{N}$ is symmetric.

For the case of yield function in the current configuration, the normal of yield surface $\bm{n}$ is

[TABLE]

$\bm{N}$ and $\bm{n}$ is related by

[TABLE]

Utilizing the transformation rule between $\bm{N}$ and $\bm{n}$, Eq.48 can be rewritten as

[TABLE]

Applying matrix exponential function, the solution of Eq.54 at step $n+1$ is

[TABLE]

By making use of Eq.44(4) $\bm{b}_{e,n}=\bm{F}\bm{C}_{p,n}^{-1}\bm{F}^{T}$, Eq.55 becomes

[TABLE]

5.2 KKT condition via variational derivation

Many physical problems involving irreversibility can be formulated as the Karush-Kuhn-Tucker (KKT) conditions [30, 11]. KKT conditions can be understood as finding the minimization of a functional subject to both equations and inequalities [11]. It generalizes the method of Lagrange multipliers by incorporating inequality constraints. In the following, we will show the KKT conditions can be derived from a functional with constraint variation.

Consider two functions $u(x)$ and $c(x)$, and we restrict that

[TABLE]

where $\epsilon$ is a small positive number. Eq.57 indicates that $\delta u$ can be varied in any direction and $\delta c$ can be done only in one direction. We may call $\delta c$ is the unidirectional variation.

A functional based on function $u(x)$ and $c(x)$ is conceptually written as

[TABLE]

The variation of $\Psi(u,c)$ reads

[TABLE]

Since we constraint that $\delta u\in[-\epsilon,\epsilon],\delta c\in[0,\epsilon]$, it is interesting to see what happens when $\delta\Psi(u,c)=0$ is enforced. $\delta\Psi(u,c)=0$ implies

[TABLE]

For the first term, $\frac{\partial\phi}{\partial u}\delta u=0,\forall\delta u$ leads to $\frac{\partial\phi}{\partial u}=0$, which is the condition of conventional variation. For the second term, the increment of $c$ is only possible when $\frac{\partial\phi}{\partial c}=0$. If $\frac{\partial\phi}{\partial c}\neq 0$, there is $\delta c=0$. Usually, the sign of $\frac{\partial\phi}{\partial c}$ can be determined in advance. If $\frac{\partial\phi}{\partial c}|_{c=0}>0$, the second equation is equivalent to

[TABLE]

which is coincidentally the same as the KKT conditions [30] if $\delta c$ is replaced with $\dot{c}$, the time derivative of $c$. In other words, the variation of an energy functional with respect to a variable of non-negative property leads to the KKT conditions.

5.3 Variational derivation of finite strain J2 plasticity and phase field model

There are several derivation of plasticity model, the maximum dissipation rule [24], the rule of dissipation inequality [57, 8], plastic metric in non-Cartesian coordinate [37]. The variational derivation presented here offers an alternative in the derivation of finite strain J2 plastic model.

Similar to the phase field model for brittle fracture, we propose the energy functional considering both the elastic energy and plastic energy and using the multiplicative elastoplasticity as

[TABLE]

where $c$ is the phase field with $c=0$ denoting no-damage material and $c=1$ denoting full-damaged material, $g(c)=(1-c)^{2}$ is the degradation function [40], $\phi^{p}(\gamma)$ is the plastic energy depending on the hardening parameter $\gamma$

[TABLE]

where $\sigma_{0}$ is the initial yield stress. In the absence of phase field, $\bm{b}_{e}=\bm{F}\bm{C}_{p}^{-1}\bm{F}^{T},\bm{\tau}=2\bm{b}_{e}\frac{\partial\phi(\bm{b}_{e})}{\partial\bm{b}_{e}}$ as shown in section.5.1. $\bm{C}_{p}^{-1}$ is the internal variable in the multiplicative elastoplasticity in section 5.1. $\phi_{+}^{e}$ and $\phi_{-}^{e}$ are the decomposition of $\phi(\bm{b}_{e})$, the stored energy energy per unit volume in the intermediate configuration, given as

[TABLE]

The variation of $\Psi$ is

[TABLE]

The second term can be simplified as

[TABLE]

Then the variation of $\Psi$ becomes

[TABLE]

Denoting \bm{\tau}(\bm{b}_{e},c)=2\big{(}g\frac{\partial\phi_{+}^{e}}{\partial\bm{b}_{e}}+\frac{\partial\phi_{-}^{e}}{\partial\bm{b}_{e}}\big{)}\bm{b}_{e} and using

[TABLE]

the last term can be simplified as

[TABLE]

Coincidentally, $||\bm{\tau}_{dev}||=\bm{\tau}(\bm{b}_{e},c):\bm{n}$ is the Von-Mises stress and $g(\sigma_{0}+h\gamma)-||\tau_{dev}||$ is equal to the yield function. It indicates that the yield function is variationally consistent.

Then the variation of $\Psi$ becomes

[TABLE]

Using integration by parts and for any $\delta c,\delta\bm{u}$, $\delta\Psi=0$ yields the governing equations

[TABLE]

Consider the fact that the $g(\sigma_{0}+h\gamma)-||\bm{\tau}_{dev}||$ is positive when the von-Mises stress is small and $\delta\gamma$ represents the plastic increment which is always non-negative, $\delta\Psi=0$ requires its last term being zero. The only possibility is that

[TABLE]

which yields Karush–Kuhn–Tucker (KKT) conditions. The derivation is similar to that in section 5.2.

Two typical decompositions of the elastic energy either based on the eigenvalue decomposition of strain tensor

[TABLE]

or the volumetric-deviatoric decomposition of the strain tensor

[TABLE]

where $\bm{\varepsilon}^{e}_{dev}=\bm{\varepsilon}^{e}-\frac{1}{3}\mbox{tr}\bm{\varepsilon}^{e}\bm{I}$.

In the current paper, considering the fact that the J2 finite strain plasticity depends on the deviatoric part of the strain tensor only, we employ the second version of the energy decomposition of strain energy density. In addition, the plasticity deformation contributes to the fracture process, the plastic work Eq.62 is added to the crack driving energy.

When large strain elastoplasticity is considered, the ratio of plastic energy is dominant compared to the elastic energy. The damage increases significantly when the plastic energy approaches the threshold value. When only the contribution of $\phi_{+}^{e}$ is considered, the fracture is brittle type, otherwise ductile type.

The damage affects both the elastic and plastic deformation. When localization occurs, undegraded stress at the crack tips increases further while stress on other parts undergo an unloading condition. So the effective Von-Mises stress depends on the damage status of that point, and the yield function is modified as

[TABLE]

where the Von-Mises stress is

[TABLE]

6 Numerical implementation

We employ the NOM scheme outlined in Section 2 to express the gradient of the displacement field. The algorithm based on spectral decomposition presented in Section 3 is used to derive and calculate the exponential matrix function, logarithmic matrix function and their first-order/second-order partial derivatives. The convergence of the numerical solution requires the local convergence of the J2 plastic model, the global convergence of the phase field and the displacement field. For the simplicity of numerical implementation, the staggered scheme in phase field model is used. The phase field model is updated using the current positive energy and historic energy as the crack driving state function. In the mechanical field, the phase field is used as a material parameter to modify the stress and the associated plastic model. Two models are calculated repeatedly until both the variations of phase field and displacement fall inside the tolerance. In the current work, we select $\epsilon_{tol}=10^{-6}$.

The J2 plastic model with phase field is achieved by using the degradation function of the von-Mises stress, see Eq.74. The actual stress is used to determine whether the yield criterion is reached. By doing so, there are two benefits: a) the actual stress located in the damaged zone is small thus avoiding the local iteration of plasticity; b) both the localizations of the damage field and plastic zone are located in a small zone. If the yield stress is degraded in the same manner, the local iteration of plasticity is required as long as the strain increases even if the material is severe damaged.

6.1 Newton-Iteration at one material point

The finite strain J2 plasticity model requires Newton iteration method to find the solution. According Section 5, the local solution of J2 plastity is equivalent to find $\bm{C}_{p}^{-1},\gamma$ at current step based on trial deformation gradient $\bm{F}$ and internal variables $\bm{C}_{p,n}^{-1},\gamma_{n}$ at previous step. The main formulas are put together as

[TABLE]

Vector $\bm{Q}$ is a nonlinear vector-valued function with unknowns specified by $\bm{h}$. The Newton-Raphson algorithm can be used to find the solution of $\bm{Q}(\bm{h})=0$ by

[TABLE]

where $\bm{h}_{k+1}$ is the vector value at iteration step $k+1$ and $\frac{\partial\bm{Q}(\bm{h}_{k})}{\partial\bm{h}_{k}},(\frac{\partial\bm{Q}(\bm{h}_{k})}{\partial\bm{h}_{k}})^{-1}$ are given by

[TABLE]

In the converged iteration of plasticity, $\bm{Q}=0$ is a function of $\bm{F}$ and $\bm{h}$. Differentiating on $\bm{Q}=0$, we have

[TABLE]

In above formula, the involved matrix expressions include $\frac{\partial\bm{h}}{\partial\bm{\bar{F}}},\frac{\partial\bm{Q}}{\partial\bar{\bm{F}}}$, which are given as follows.

[TABLE]

In sum, the algorithm to implement the local iteration and consistent tangent stiffness matrix are summarized in Algorithm 1.

Algorithm 1: Local iteration for J2 plasticity

$\displaystyle\mathcal{F}^{trial}=\mathcal{F}(\bm{F},\bm{h}_{n})$

$\displaystyle\mbox{IF }\mathcal{F}^{trial}<0,\mbox{ then in elastic state}$

$\displaystyle\quad\bm{h}=\bm{h}_{n}$

ELSE: in plastic state, begin local plastic iteration

$\displaystyle\quad\bm{h}=\bm{h}_{n}$

REPEAT

$\displaystyle\quad\quad\bm{A}=\frac{\partial\bm{Q}(\bm{F},\bm{h},\bm{h}_{n})}{\partial\bm{h}}$

$\displaystyle\quad\quad\Delta\bm{h}=-\bm{A}^{-1}\bm{Q}$

$\displaystyle\quad\quad\bm{h}=\bm{h}+\Delta\bm{h}$

$\displaystyle\quad\mbox{ UNTIL }\|\Delta\bm{h}\|\leq\varepsilon_{h}$

$\displaystyle\quad\bm{h}=\bm{h},\frac{\partial\bm{h}}{\partial\bm{F}}=-\bm{A}(\bm{F},\bm{h},\bm{h}_{,n})^{-1}\frac{\partial\bm{Q}(\bm{F},\bm{h},\bm{h}_{,n})}{\partial\bm{F}}$

END IF

Let $\bm{F}$ be written in vector form $\bar{\bm{F}}$ as

[TABLE]

Let the list of unknown displacement in support denoted by $\bm{U}$, $\bar{\bm{F}}$ is a function of $\bm{U}$. According to Eq.14, we have $\frac{\partial\bar{\bm{F}}}{\partial\bm{U}}=\bm{B}$.

Then the local residual vector and consistent tangent stiffness matrix of one particle based on unknown vector $\bm{U}$ are calculated by

[TABLE]

[TABLE]

where $\Delta V$ is the volume of the material point. After assembling $\mathbb{R}_{i}$ and $\mathbb{K}_{i}$ of all particles into the global residual and global stiffness matrix, the standard Newton-Iteration algorithm can be applied.

The numerical method in the paper is implemented using the nonlocal operator method. The number of particles in support is chosen as 9 in 2D and 27 in 3D, where the support size is determined by the furthest particles in support for each particle. The weight function is selected as $w(r)=1/r^{2}$. The nonlocal gradient in discrete form plays the same role as the shape function in finite element method. Since nodal integral is used, the hourglass energy functional is added to the weak form to suppress the zero-energy mode. The calculation of residual and tangent stiffness matrix for hourglass energy function is based on Eq.17.

7 Numerical examples

7.1 Phase field based on Hencky strain without plasticity: Single-edge-notched tension test

The current phase field scheme can recover the small strain phase field model when the $G_{c}$ in small strain brittle fracture is used. The plasticity deformation can be prevented by making the yield stress large enough. In this subsection, we model the single-edge-notched tension test, which is a squared plate with initial notched crack as shown in Fig.2. The material parameters are set as $\lambda=121.1538$ kN/mm2 and $\mu=80.7692$ kN/mm2 for elastic constants, $G_{c}=2.7\times 10^{-3}$ kN/mm for critical energy release rate. These parameters are identical to that used in the small strain brittle fracture phase field in Ref [43]. Two displacement conditions are tested: Case a) for tensile boundary condition and Case b) for shear boundary conditions. The plate is discretized into $100\times 100$ material points, the phase field length scale is selected as $l=10^{-2}$ mm. The displacement load is monotonic applied with fixed displacement increment $\Delta u=1\times 10^{-5}$ mm.

The evolution of the phase field for tensile tests with a discretization of $200\times 200$ is depicted in Fig.4. The final displacement field in y-direction is given in Fig.4. The displacement-load curve is plotted in Fig.5.

The phase field evolution of shear test for brittle fracture is depicted in Fig.6. The displacement field is given in Fig.7. When plasticity is not involved and the critical energy release rate for brittle fracture is used, the Hencky strain phase field model can recover the small strain version. The load curve of notched plate modeled by the current scheme is compared to the traditional small strain finite element method in Fig.8. The phase field length scale dependence is observed in the current simulation. On the other hand, the Hencky strain formulation considers the geometric nonlinearity and shows a slightly smooth degradation.

7.2 Phase field based on Hencky strain elastoplasticity: Single-edge-notched test with ductile fracture

The plate is discretized into 200x200 particles. The ductile critical energy release rate is selected as $G_{c}^{d}=(5G_{c},10G_{c},15G_{c})$, where $G_{c}=2.7\times 10^{-3}$ kN/mm is the critical energy release for brittle fracture. The Hencky strain J2 plasticity is used. For the case of $G_{c}^{d}=5G_{c}$, the final phase field, displacement field, equivalent plasticity and deformed configuration are provided in Fig.9. It can be observed that the crack propagates along the direction of the initial crack and the plastic deformation and ductile fracture happens in the same region. The plastic energy drives the propagation of the phase field dominantly. Conversely, the evolution of the phase field unloads the internal force in other regions, which induces a strain localization and large plastic deformation only in the damaged region. Compared with the brittle phase field, the notched plate undergoes a large deformation before the phase field becomes significant. The reaction forces for different $G_{c}^{d}$ is provided in Fig.10. The displacement-load decreases gradually due to the dissipation in plastic deformation.

The evolution of phase field for different time instants is illustrated in Fig.11 and the equivalent plasticity field is shown in Fig.12. The tensile load results in a large deformation along the path of phase field region. It is observed that the large plastic region and the phase field region coincide, which shows that the energy due to plastic deformation contributes significantly to the development of phase field.

7.3 Isotropic hardening J2 plasticity: necking of a rod

Necking of a rod is a classical benchmark problem for isotropic elastic-plastic material subjected to finite plasticity, see for example [57, 47, 38]. The length of the rod in the initial configuration is $l=53.34$ mm, the radius $r_{0}=6.4135$ mm. A finer discretization is employed in the middle of the rod close to the necked zone. In order to trigger a necking, an imperfection of the center of the rod is introduced by reducing the radius to $r=0.982r_{0}$ gradually. The nonlinear isotropic hardening response is described by the yielding function

[TABLE]

whose material parameters are given in Table 1.

The number of nodes is 4758 ($\approx 5k$) for coarse mesh as shown in Fig.13 and 49320 ($\approx 50k$) for refined mesh.

The load-displacement is plotted and compared to the finite element method, as shown in Fig.15. The results depict that NOM is suitable for large strain J2 plasticity models and can achieve similar results as finite element methods.

7.4 Necking of a rod with ductile fracture

The material parameters, discretization of the 3D rod and the boundary conditions are the same as those in Section 7.3. In order to consider the evolution of ductile fracture, the phase field model is activated. The phase field length scale is selected as the average particle size of $l=1.2$ mm and the critical energy release rate as $G_{c}=2.7\times 10^{-3}$ kN/mm. We study the resultant reaction forces under the influence of critical release energy rate. For the case of $25G_{c}$, the actually deformed particle distribution is given in Fig.17. Although the NOM scheme used here is based on nodal discretization, an auxiliary mesh is employed to show the deformation field and other physical quantities such as plasticity and phase field. The deformed configuration, phase field distribution and equivalent plasticity is shown in Fig.16.

The displacement-load curve for ductile fracture of a rod for different critical energy release rate is depicted in Fig.18. Conventional J2 plasticity is not easy to describe the damaging process. With the aid of the phase field method, the ductile fracture considering the contribution of phase field deformation can be modeled. It is observed that the ductile fracture differs from the brittle fracture by undergoing firstly a plastic deformation. Since the elastic energy occupies a small part of the total energy, the accumulated plastic energy reaches the threshold value, the damage occurs immediately and the bearing capacity of the rod reduces very fast. With the increase of the critical energy release rate, the critical displacement increases proportionally.

8 Conclusions

In this paper, we have developed a fintie strain plastic model with phase field method for the modeling of ductile fractures. The model contains several parts: the J2 plasticity model, the consistent tangent stiffness matrix based on matrix functions, the variational derivation of phase field plastic model based on its energy form. The implementation is based on nonlocal operator method, which uses the nonlocal derivative to replace the gradient of shape functions.

Several numerical examples including the brittle fracture subjected to tensile load and shear load, ductile fracture in 2D and fintie strain J2 plasticity in 3D, and ductile fracture in 3D are presented to show the capability of current scheme.

Acknowledgements

The first author acknowledges the financial support of the EU project under the title of “Computational Modeling, Topological Optimization and Design of Flexoelectric Nano Energy Harvesters” (ERC COTOFLEXI 802205).