Thermomechanical Peridynamic Modeling for Ductile Fracture

1. Introduction

Many ductile materials, such as superalloys, are subjected to high temperatures for long periods, resulting in irreversible plastic deformation and even ductile fracture. Although researchers have conducted a number of numerical studies on the failure of ductile materials [1,2,3], the methods proposed are not suitable for predicting the process of crack initiation and propagation. A reason for this is that the constitutive models used and the corresponding numerical simulations were based on classical continuum mechanics (CCM). When CCM is used for fracture simulations, one inevitably encounters the singularity problem posed by discontinuities. To avoid the difficulties repeatedly encountered by CCM, a nonlocal solid mechanics theory called peridynamics (PD) was proposed by Silling in 2000 [4]. In recent years, the advantages of metal fracture modeling and numerical simulation based on PD are gradually becoming apparent [5,6].

PD assumes that there still are interactions between non-contact material points separated by a finite distance, and describes the mechanical behaviors of materials by solving integral equations of nonlocal effects in space, avoiding the singularity that occurs when CCM faces the fracture problem [7,8]. Therefore, PD models are suitable for simulating the fracture process and do not require complex crack propagation criteria [9].

However, nonlocal integration greatly increases calculation costs, and traditional traction boundary conditions cannot be applied directly in PD models. An effective way to solve this problem is to couple PD with CCM [10,11,12,13,14]. Lubineau et al. proposed a Morphing method to couple PD with CCM based on the principle of energy density balance [15], Azdoud et al. used this method to study quasi-static fracture problems [16]. Han et al. further developed a new model of coupling damage mechanics with PD through the Morphing method and used it to simulate the entire failure process of structures from elastic deformation and damage to complete fracture [17].

PD can be used to simulate not only the behavior of linear elastic materials but also the plastic deformation of materials. Macek and Silling proposed a perfectly elastoplastic bond constitutive model by permitting permanent deformation of individual bonds [18]. Zhan et al. proposed a rate-dependent plastic PD model to simulate the dynamic response of particle-reinforced metal matrix composites [19]. Madenci and Oterkus proposed an ordinary state-based PD model for plastic deformation based on von Mises yield criteria [20]. Liu et al. proposed an ordinary state-based PD model for nonlinear hardening plastic materials [21]. Although great progress has been made in the PD modeling of plastic materials, the above models did not take temperature into account. Both thermal and plastic deformation should be considered in the simulation of the plastic fracture of certain materials, such as superalloys.

In this paper, a plastic bond constitutive model is developed based on a thermoelastic bond-based PD model to describe the deformation and fracture of superalloys at high temperatures. The Morphing method is used to establish a thermoelastic model coupling PD with CCM to limit the plastic PD model to the fracture region, thereby reducing the computational cost.

The remainder of this paper is organized as follows: Section 2 proposes a new bond-based PD constitutive model by considering plasticity and thermal expansion in bond deformation. Section 3 derives a thermoelastic coupling model between PD and CCM through the Morphing method. Then, in Section 4, we describe an iterative algorithm of plastic-fracture calculations. We present numerical results for the fracture process in a superalloy structure and compare them with experimental observations in Section 5. Finally, Section 6 concludes the paper.

2. Thermal Elastoplastic Bond-based Peridynamic Model

PD theory assumes that there is an interaction between non-contact material points separated by a finite distance. Therefore, in the reference configuration $\mathsf{\Omega }$ of the object, the interaction between any point $\mathit{x}$ and a point ${\mathit{x}}^{\prime }$ in its neighborhood is described by the pairwise force function $\mathit{f}$, and the balance equation at point $\mathit{x}$ can be written as

$${\int }_{{\mathcal{H}}_{\delta }\left(\mathit{x}\right)}\mathit{f}({\mathit{x}}^{\prime },\mathit{x})d{V}_{{\mathit{x}}^{\prime }}+\mathit{b}\left(\mathit{x}\right)=\mathbf{0}$$

(1)

where $\mathit{b}$ is a prescribed body force density and ${\mathcal{H}}_{\delta }\left(\mathit{x}\right)\subset \mathsf{\Omega }$ is a circular region with $\mathit{x}$ as its center and $\delta$ as its radius (see Figure 1). Thus, $\delta$ is referred to as the horizon [4].

The pairwise force function $\mathit{f}$ is defined as

$$\mathit{f}({\mathit{x}}^{\prime },\mathit{x})=\hat{\mathit{f}}({\mathit{x}}^{\prime },\mathit{x})-\widehat{\mathit{f}}(\mathit{x},{\mathit{x}}^{\prime })$$

(2)

where $\widehat{\mathit{f}}({\mathit{x}}^{\prime },\mathit{x})$ is the partial interaction due to the action of point ${\mathit{x}}^{\prime }$ over point $\mathit{x}$ and, correspondingly, $\widehat{\mathit{f}}(\mathit{x},{\mathit{x}}^{\prime })$ is the partial interaction due to the action of point $\mathit{x}$ over point ${\mathit{x}}^{\prime }$. By introducing plastic deformation and thermal expansion into the bond deformation calculation, the partial interaction can be written as

$$\widehat{\mathit{f}}({\mathit{x}}^{\prime },\mathit{x})=\frac{1}{2}c(\mathit{x},\mathit{\xi })\left({\eta }_{\mathit{\xi }}({\mathit{x}}^{\prime },\mathit{x})-{\overline{\eta }}_{\mathit{\xi }}-a(\mathit{x},\mathit{\xi })\widehat{T}(\mathit{x},\mathit{\xi })\right){\mathit{e}}_{\mathit{\xi }}$$

(3)

where $\mathit{\xi }={\mathit{x}}^{\prime }-\mathit{x}$ is the relative position vector herein referred to as “bond", ${\eta }_{\mathit{\xi }}({\mathit{x}}^{\prime },\mathit{x})=\mathit{\eta }({\mathit{x}}^{\prime },\mathit{x})·{\mathit{e}}_{\mathit{\xi }}$ denotes the projection of the bond deformation at point $\mathit{x}$ over the bond, $\mathit{\eta }({\mathit{x}}^{\prime },\mathit{x})=\mathit{u}\left({\mathit{x}}^{\prime }\right)-\mathit{u}\left(\mathit{x}\right)$, and ${\mathit{e}}_{\mathit{\xi }}=\raisebox{1ex}{$\mathit{\xi }$}\!\left/ \!\raisebox{-1ex}{$\left|\mathit{\xi }\right|$}\right.$. In addition, $\mathit{u}$ is the displacement and ${\mathit{e}}_{\mathit{\xi }}$ is the unit direction vector of the bond. In Eq.(3), $a(\mathit{x},\mathit{\xi })\widehat{T}(\mathit{x},\mathit{\xi })$ represents the deformation of the bond due to thermal expansion [22], where $a(\mathit{x},\mathit{\xi })$ is the coefficient of thermal expansion of the bond in the PD model. We define $\widehat{T}(\mathit{x},\mathit{\xi })$ as the effective temperature of the bond. In other words, $\widehat{T}$ is determined by $T\left(\mathit{x}\right)$ and $T\left({\mathit{x}}^{\prime }\right)$, which are the temperature variations between the current temperatures and the reference temperatures of points $\mathit{x}$ and ${\mathit{x}}^{\prime }$, respectively. Thus, we have

$$\widehat{T}(\mathit{x},\mathit{\xi })=\widehat{T}({\mathit{x}}^{\prime },-\mathit{\xi })$$

(4)

Moreover, ${\overline{\eta }}_{\mathit{\xi }}=\overline{\mathit{\eta }}·{\mathit{e}}_{\mathit{\xi }}$, where $\overline{\mathit{\eta }}$ is the historical plastic elongation of the bond defined by:

$$\overline{\mathit{\eta }}\left(0\right)=0,\dot{\overline{\mathit{\eta }}}=\left\{\begin{array}{cc}\dot{\mathit{\eta }}\hfill & \mathrm{if}\frac{\left|\mathit{\xi }+\mathit{\eta }-\overline{\mathit{\eta }}\right|-\left|\mathit{\xi }\right|}{\left|\mathit{\xi }\right|}\ge s_Y\hfill \\ 0\hfill & \mathrm{otherwise}\hfill \end{array}\right.$$

(5)

where $s_Y$ is the yield stretch. It can be seen that, at the bond level, this material is elastic-perfectly plastic. However, for the entire structure, the material model has a strain hardening effect because not all bonds will yield at a certain moment or under deformation. It is worth noting that $s_Y$ is an intrinsic parameter of the material, and its value is related to its engineering ultimate strength [18]. Moreover, the fracture behavior of materials can be described via bond failure [23]. It is assumed that the bond will break irreversibly if its stretch s exceeds the critical value $s_C$ (generally, $s_C>s_Y$), where $s=\frac{\left|\mathit{\xi }+\mathit{\eta }\right|-\left|\mathit{\xi }\right|}{\left|\mathit{\xi }\right|}$.

After introducing a plastic elongation at the bond level, we then give an average measure of non-local plastic deformation at a point, which can be expressed as:

$${\varphi }_p\left(\mathit{x}\right)=\frac{{\int }_{{\mathcal{H}}_{\delta }\left(\mathit{x}\right)}\overline{s}d{V}_{\mathit{\xi }}}{{\int }_{{\mathcal{H}}_{\delta }\left(\mathit{x}\right)}d{V}_{\mathit{\xi }}}$$

(6)

where $\overline{s}=\frac{|\mathit{\xi }+\overline{\eta }|-|\mathit{\xi }|}{\left|\mathit{\xi }\right|}$ is the plastic stretch of a bond.

3. Thermoelastic Hybrid Peridynamics and Classical Continuum Mechanics (PD-CCM) Model

Hybrid PD-CCM models retain the advantages of PD in discontinuities and improve calculation efficiency [10,11,12,13,14]. In this section, we adopt the Morphing method [15] to establish a thermoelastic hybrid PD-CCM model.

3.1. Governing Equations

We divide the whole solution region into three non-overlapping subregions, i.e., $\mathsf{\Omega }={\mathsf{\Omega }}_1\cup {\mathsf{\Omega }}_2\cup {\mathsf{\Omega }}_m$, and ${\mathsf{\Omega }}_1\cap {\mathsf{\Omega }}_2=\varnothing$, ${\mathsf{\Omega }}_1\cap {\mathsf{\Omega }}_m=\varnothing$, ${\mathsf{\Omega }}_2\cap {\mathsf{\Omega }}_m=\varnothing$. The PD model is confined to subregion ${\mathsf{\Omega }}_2$, where plastic deformation and fracture may occur, and the CCM model is employed in subregions ${\mathsf{\Omega }}_1$ to reduce computational costs. The two models transition through subregion ${\mathsf{\Omega }}_m$. As shown in Figure 2, the displacement boundary condition $\overline{\mathit{u}}$ is applied to part ${\mathsf{\Gamma }}_{\mathit{u}}$ of $\partial \mathsf{\Omega }$, and the traction load $\overline{\mathit{F}}$ is applied to part ${\mathsf{\Gamma }}_{\mathit{F}}$ of $\partial \mathsf{\Omega }$. $\mathit{n}$ is a unit vector representing the outer normal, and the entire solution region $\mathsf{\Omega }$ is affected by the body force $\mathit{b}$ and the steady-state temperature field $T\left(\mathit{x}\right)$.

The Morphing method is based on a set of governing equations that are valid for any point in $\mathsf{\Omega }$, and it defines the material properties of the two models as functions that change with position [24]. The governing equations are as follows:

• Kinematic admissibility and compatibility

  $$\begin{array}{cccc}\hfill \mathit{\epsilon }\left(\mathit{x}\right)& =\frac{1}{2}\left(\nabla \mathit{u}\left(\mathit{x}\right)+{(\nabla \mathit{u})}^T\left(\mathit{x}\right)\right)\hfill & \hfill & \forall \mathit{x}\in \mathsf{\Omega }\setminus {\mathsf{\Omega }}_2\hfill \end{array}$$

  (7)

  $$\begin{array}{cccc}\hfill {\eta }_{\mathit{\xi }}\left({\mathit{x}}^{\prime },\mathit{x}\right)& =u_{\mathit{\xi }}\left({\mathit{x}}^{\prime }\right)-u_{\mathit{\xi }}\left(\mathit{x}\right)\hfill & \hfill & \forall \mathit{x},{\mathit{x}}^{\prime }\in \mathsf{\Omega }\setminus {\mathsf{\Omega }}_1\hfill \end{array}$$

  (8)

  $$\begin{array}{cccc}\hfill \mathit{u}\left(\mathit{x}\right)& =\overline{\mathit{u}}\left(\mathit{x}\right)\hfill & \hfill & \forall \mathit{x}\in {\mathsf{\Gamma }}_{\mathit{u}}\hfill \end{array}$$

  (9)
• Static admissibility

  $$\begin{array}{cccc}\hfill div\mathit{\sigma }\left(\mathit{x}\right)+{\int }_{{\mathcal{H}}_{\delta }\left(\mathit{x}\right)}\mathit{f}({\mathit{x}}^{\prime },\mathit{x})d{V}_{{\mathit{x}}^{\prime }}& =-\mathit{b}\left(\mathit{x}\right)\hfill & \hfill & \forall \mathit{x}\in \mathsf{\Omega }\hfill \end{array}$$

  (10)

  $$\begin{array}{cccc}\hfill \mathit{\sigma }\left(\mathit{x}\right)·\mathit{n}\left(\mathit{x}\right)& =\overline{\mathit{F}}\hfill & \hfill & \forall \mathit{x}\in {\mathsf{\Gamma }}_{\mathit{F}}\hfill \end{array}$$

  (11)
• Constitutive equations

  $$\begin{array}{cccc}\hfill \mathit{\sigma }\left(\mathit{x}\right)& =\mathit{E}\left(\mathit{x}\right):\mathit{\epsilon }\left(\mathit{x}\right)-\mathit{\beta }\left(\mathit{x}\right)T\left(\mathit{x}\right)\hfill & \hfill & \forall \mathit{x}\in \mathsf{\Omega }\setminus {\mathsf{\Omega }}_2\hfill \\ \hfill \mathit{f}\left({\mathit{x}}^{\prime },\mathit{x}\right)& =\frac{1}{2}\left(c(\mathit{x},\mathit{\xi })+c\left({\mathit{x}}^{\prime },\mathit{\xi }\right)\right)\left({\eta }_{\mathit{\xi }}\left({\mathit{x}}^{\prime },\mathit{x}\right)-{\overline{\eta }}_{\mathit{\xi }}\right){\mathit{e}}_{\mathit{\xi }}\hfill & \hfill & \hfill \end{array}$$

  (12)

  $$\begin{array}{cccc}\hfill & -\frac{1}{2}\left(b(\mathit{x},\mathit{\xi })\widehat{T}(\mathit{x},\mathit{\xi })+b\left({\mathit{x}}^{\prime },\mathit{\xi }\right)\widehat{T}({\mathit{x}}^{\prime },-\mathit{\xi })\right){\mathit{e}}_{\mathit{\xi }}\hfill & \hfill & \forall \mathit{x}\in \mathsf{\Omega }\setminus {\mathsf{\Omega }}_1\hfill \end{array}$$

  (13)

Here, $\mathit{\sigma }$ and $\mathit{\epsilon }$ are the stress and strain tensors, respectively. In Eq.(12), $\mathit{E}\left(\mathit{x}\right)$ is the stiffness tensor, and $\mathit{\beta }\left(\mathit{x}\right)=\mathit{E}\left(\mathit{x}\right):\mathit{\alpha }\left(\mathit{x}\right)$ is the thermal modulus, where $\mathit{\alpha }\left(\mathit{x}\right)$ is the thermal expansion coefficient tensor.

For homogeneous materials, the PD micromodulus function $c(\mathit{x},\mathit{\xi })$ and the thermal expansion coefficient function $a(\mathit{x},\mathit{\xi })$ can be redefined as:

$$c(\mathit{x},\mathit{\xi })=w\left(\mathit{x}\right)c^0\left(\right|\mathit{\xi }\left|\right)$$

(14)

$$a(\mathit{x},\mathit{\xi })=v\left(\mathit{x}\right)a^0\left(\right|\mathit{\xi }\left|\right)$$

(15)

where $w\left(\mathit{x}\right)$ and $v\left(\mathit{x}\right)$ are the Morphing functions, $0\le w\left(\mathit{x}\right),v\left(\mathit{x}\right)\le 1,\forall \mathit{x}\in \mathsf{\Omega }$. $c^0\left(\right|\mathit{\xi }\left|\right)$ and $a^0\left(\right|\mathit{\xi }\left|\right)$ are the PD micromodulus and thermal expansion coefficient, respectively. We define the PD thermal modulus as $b^0\left(\right|\mathit{\xi }\left|\right)=c^0\left(\right|\mathit{\xi }\left|\right)a^0\left(\right|\mathit{\xi }\left|\right)$. Then, the pairwise force function $\mathit{f}$ can be rewritten as in Eq.(13), where

$$b(\mathit{x},\mathit{\xi })=w\left(\mathit{x}\right)v\left(\mathit{x}\right)b^0\left(\right|\mathit{\xi }\left|\right)$$

(16)

3.2. Stiffness/Thermal Modulus Constraint Equations

The overall governing equations of the hybrid model were given in Section 3.1. It can be seen that the functions $w\left(\mathit{x}\right)$, $v\left(\mathit{x}\right)$ and $\mathit{E}\left(\mathit{x}\right)$, $\mathit{\beta }\left(\mathit{x}\right)$ jointly guide the transition between the CCM and PD model. $w\left(\mathit{x}\right)$ and $v\left(\mathit{x}\right)$ are the functions that determine the division of the hybrid model, so how we define $\mathit{E}\left(\mathit{x}\right)$ and $\mathit{\beta }\left(\mathit{x}\right)$ is key to successfully couple both models together. In this subsection, we derive the constraint functions of $\mathit{E}\left(\mathit{x}\right)$ and $\mathit{\beta }\left(\mathit{x}\right)$ based on the point-by-point equivalency of energy density.

The plastic deformation is confined to the pure peridynamic region, that is, ${\overline{\eta }}_{\mathit{\xi }}\equiv 0$ outside this region. Thus, the constitutive equation in Equation (13) can be simplified to an elastic version. Let us consider a subregion ${\mathsf{\Omega }}_0=\{\mathit{x}\in \mathsf{\Omega }:{\mathcal{H}}_{\delta }\left(\mathit{x}\right)\subset \mathsf{\Omega }\}$, such that we still have a complete integral region when the PD model is used near the boundary. We can derive the following equations according to previous work [24]. The following assertions hold true for any point $\mathit{x}\in {\mathsf{\Omega }}_0$:

• If and only if $\mathit{E}\left(\mathit{x}\right)={\mathit{E}}^0$, $\mathit{\beta }\left(\mathit{x}\right)={\mathit{\beta }}^0$ and $w\left({\mathit{x}}^{\prime }\right),v\left({\mathit{x}}^{\prime }\right)=0$, $\forall {\mathit{x}}^{\prime }\in {\mathcal{H}}_{\delta }\left(\mathit{x}\right)$, the model is restricted to the pure CCM model at point $\mathit{x}$. Then the elastic energy density at point $\mathit{x}$ can be written as

  $$W\left(\mathit{x}\right)=\frac{1}{2}\mathit{\epsilon }\left(\mathit{x}\right):{\mathit{E}}^0:\mathit{\epsilon }\left(\mathit{x}\right)-\mathit{\epsilon }\left(\mathit{x}\right):{\mathit{\beta }}^{0}T\left(\mathit{x}\right)+W^0\left(T\right)$$

  (17)
• If and only if $\mathit{E}\left(\mathit{x}\right)=0$, $\mathit{\beta }\left(\mathit{x}\right)=0$ and $w\left({\mathit{x}}^{\prime }\right),v\left({\mathit{x}}^{\prime }\right)=1$, $\forall {\mathit{x}}^{\prime }\in {\mathcal{H}}_{\delta }\left(\mathit{x}\right)$, the model is restricted to the pure PD model at point $\mathit{x}$. Considering Equation (4), the elastic energy density at point $\mathit{x}$ can be written as

  $$W\left(\mathit{x}\right)=\frac{1}{4}{\int }_{{\mathcal{H}}_{\delta }\left(\mathit{x}\right)}{c}^0\left(\right|\mathit{\xi }\left|\right){\eta }_{\mathit{\xi }}^2({\mathit{x}}^{\prime },\mathit{x})d{V}_{{\mathit{x}}^{\prime }}-\frac{1}{2}{\int }_{{\mathcal{H}}_{\delta }\left(\mathit{x}\right)}{b}^0\left(\right|\mathit{\xi }\left|\right)\widehat{T}(\mathit{x},\mathit{\xi }){\eta }_{\mathit{\xi }}({\mathit{x}}^{\prime },\mathit{x})d{V}_{{\mathit{x}}^{\prime }}+W^0\left(T\right)$$

  (18)
• If and only if $\mathit{E}\left(\mathit{x}\right)\ne 0$, $\mathit{\beta }\left(\mathit{x}\right)\ne 0$ and $\exists {\mathit{x}}^{\prime }\in {\mathcal{H}}_{\delta }\left(\mathit{x}\right)$ such that $0<w\left(\mathit{x}\right)<1,0<v\left(\mathit{x}\right)\le 1$, the hybrid model must be used at point $\mathit{x}$. Considering Equation (4), the elastic energy density at point $\mathit{x}$ can be written as

  $$\begin{array}{c}\hfill W\left(\mathit{x}\right)=\frac{1}{2}\mathit{\epsilon }\left(\mathit{x}\right):\mathit{E}\left(\mathit{x}\right):\mathit{\epsilon }\left(\mathit{x}\right)+\frac{1}{4}{\int }_{{\mathcal{H}}_{\delta }\left(\mathit{x}\right)}{c}^0\left(\right|\mathit{\xi }\left|\right)\frac{w\left({\mathit{x}}^{\prime }\right)+w\left(\mathit{x}\right)}{2}{\eta }_{\mathit{\xi }}^2({\mathit{x}}^{\prime },\mathit{x})d{V}_{{\mathit{x}}^{\prime }}+W^0\left(T\right)\\ \hfill -\mathit{\epsilon }\left(\mathit{x}\right):\mathit{\beta }\left(\mathit{x}\right)T\left(\mathit{x}\right)-\frac{1}{2}{\int }_{{\mathcal{H}}_{\delta }\left(\mathit{x}\right)}{b}^0\left(\right|\mathit{\xi }\left|\right)\frac{w\left({\mathit{x}}^{\prime }\right)v\left({\mathit{x}}^{\prime }\right)+w\left(\mathit{x}\right)v\left(\mathit{x}\right)}{2}\widehat{T}(\mathit{x},\mathit{\xi }){\eta }_{\mathit{\xi }}({\mathit{x}}^{\prime },\mathit{x})d{V}_{{\mathit{x}}^{\prime }}\end{array}$$

  (19)

Remark: 

$W^0\left(T\right)$ corresponds to the $T^2$ term in the elastic energy density equations.

For homogeneous materials, the elastic energy density at any point should be independent of the definition of the Morphing functions $w\left(\mathit{x}\right)$ and $v\left(\mathit{x}\right)$ and only depend on the displacement and temperature field. This means that the elastic energy densities calculated at the same point with different expressions are equivalent to each other, and thus the corresponding terms of force and temperature must be equal. Therefore, the terms of Equation (17) are equal to the corresponding ones in Equation (18) respectively, i.e.,

$$\begin{array}{cc}\hfill \frac{1}{2}\mathit{\epsilon }\left(\mathit{x}\right):{\mathit{E}}^0:\mathit{\epsilon }\left(\mathit{x}\right)& =\frac{1}{4}{\int }_{{\mathcal{H}}_{\delta }\left(\mathit{x}\right)}{c}^0\left(\right|\mathit{\xi }\left|\right){\eta }_{\mathit{\xi }}^2({\mathit{x}}^{\prime },\mathit{x})d{V}_{{\mathit{x}}^{\prime }}\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill \mathit{\epsilon }\left(\mathit{x}\right):{\mathit{\beta }}^{0}T\left(\mathit{x}\right)& =\frac{1}{2}{\int }_{{\mathcal{H}}_{\delta }\left(\mathit{x}\right)}{b}^0\left(\right|\mathit{\xi }\left|\right)\widehat{T}(\mathit{x},\mathit{\xi }){\eta }_{\mathit{\xi }}({\mathit{x}}^{\prime },\mathit{x})d{V}_{{\mathit{x}}^{\prime }}\hfill \end{array}$$

(21)

Furthermore, we assumed that both temperature and elastic deformation in the small neighborhood of point $\mathit{x}$ are uniform such that

$$\mathit{\epsilon }\left({\mathit{x}}^{\prime }\right)\simeq \mathit{\epsilon }\left(\mathit{x}\right)=\overline{\mathit{\epsilon }}\mathrm{and}{\eta }_{\mathit{\xi }}\left({\mathit{x}}^{\prime },\mathit{x}\right)=\frac{\mathit{\xi }·\overline{\mathit{\epsilon }}·\mathit{\xi }}{\left|\mathit{\xi }\right|},\forall {\mathit{x}}^{\prime }\in {\mathcal{H}}_{\delta }\left(\mathit{x}\right)$$

(22)

$$\widehat{T}(\mathit{x},\mathit{\xi })\simeq T\left(\mathit{x}\right)=\overline{T},\forall {\mathit{x}}^{\prime }\in {\mathcal{H}}_{\delta }\left(\mathit{x}\right)$$

(23)

Thus, Equations (20) and (21) leads to:

$$\begin{array}{cc}\hfill {\mathit{E}}^0& =\frac{1}{2}{\int }_{{\mathcal{H}}_{\delta }\left(\mathit{x}\right)}{c}^0\left(\right|\mathit{\xi }\left|\right)\frac{\mathit{\xi }\otimes \mathit{\xi }\otimes \mathit{\xi }\otimes \mathit{\xi }}{{\left|\mathit{\xi }\right|}^2}dV{\mathit{x}}^{\prime }\hfill \end{array}$$

(24)

$$\begin{array}{cc}\hfill {\mathit{\beta }}^0& =\frac{1}{2}{\int }_{{\mathcal{H}}_{\delta }\left(\mathit{x}\right)}{b}^0\left(\right|\mathit{\xi }\left|\right)\frac{\mathit{\xi }\otimes \mathit{\xi }}{\left|\mathit{\xi }\right|}dV{\mathit{x}}^{\prime }\hfill \end{array}$$

(25)

Similarly, the terms of Equation (17) are equal to the corresponding ones in Equation (19), and by considering Equations (22)–(25) we get

$$\begin{array}{cc}\hfill \mathit{E}\left(\mathit{x}\right)& ={\mathit{E}}^0-\frac{1}{2}{\int }_{{\mathcal{H}}_{\delta }\left(\mathit{x}\right)}{c}^0\left(\right|\mathit{\xi }\left|\right)\frac{w\left({\mathit{x}}^{\prime }\right)+w\left(\mathit{x}\right)}{2}\frac{\mathit{\xi }\otimes \mathit{\xi }\otimes \mathit{\xi }\otimes \mathit{\xi }}{{\left|\mathit{\xi }\right|}^2}dV{\mathit{x}}^{\prime }\hfill \end{array}$$

(26)

$$\begin{array}{cc}\hfill \mathit{\beta }\left(\mathit{x}\right)& ={\mathit{\beta }}^0-\frac{1}{2}{\int }_{{\mathcal{H}}_{\delta }\left(\mathit{x}\right)}{b}^0\left(\right|\mathit{\xi }\left|\right)\frac{w\left({\mathit{x}}^{\prime }\right)v\left({\mathit{x}}^{\prime }\right)+w\left(\mathit{x}\right)v\left(\mathit{x}\right)}{2}\frac{\mathit{\xi }\otimes \mathit{\xi }}{\left|\mathit{\xi }\right|}dV{\mathit{x}}^{\prime }\hfill \end{array}$$

(27)

In summary, the thermoelastic hybrid PD-CCM model is established by Equations (7)–(16), (26) and (27). The partitioning of the hybrid model is determined by the Morphing functions $w\left(\mathit{x}\right)$ and $v\left(\mathit{x}\right)$, and the definitions of $w\left(\mathit{x}\right)$ and $v\left(\mathit{x}\right)$ are independent of each other.

4. Plastic-Fracture Calculations

In this paper, the constraint equations of the hybrid model are used to ensure that PD applies only in the critical region where plastic deformation and cracks may occur. The plastic bond constitutive equation defined in Equation (3) is used to perform iterative calculations.

Figure 3 shows the flowchart of the implicit iterative algorithm used for plastic-fracture simulation. In this algorithm, the coupling zone and the CCM zone are discretized by continuous finite elements, whereas the pure PD zone is discretized by discrete finite elements [25]. Discrete elements allow cracks to grow freely along their boundaries. The boundary conditions are imposed gradually in N increments. In each incremental step, the elongation of each bond is calculated according to the displacement field obtained by solving the corresponding linear equations. If any bond yields, the residual plastic force needs to be calculated to update the force vector at the right side of the equations. This step is repeated until the iterative solution of the system of linear equations converges.

After the plastic iterative calculation converges, if any bonds break at the current step, it is necessary to update the stiffness matrix by removing the contribution of these bonds. Change in the stiffness matrix will inevitably alter the calculated displacement field, so it is necessary to recalculate the residual plastic force and solve the equations until the plastic iteration converges and no bonds break anymore at the current step. Afterwards, the step is incremented and calculations are repeated until the last load step is completed.

5. Numerical Examples

This section presents two numerical examples to validate the proposed model. The first is a two-dimensional simplified model of an actual experiment involving a GH4099 superalloy structure under high-temperature environments. We compare the simulated fracture modes and force-displacement curves with the experimental results to assess the ability of the proposed model for simulating the failure of ductile materials. The second numerical example involves a four-point bending beam under uniform and non-uniform temperature fields to investigate the effects of temperature field and plastic deformation on crack propagation.

5.1. Example 1: GH4099 Superalloy Structure

The geometry of the experimental GH4099 superalloy specimen is shown in Figure 4. As shown in Figure 4b, the top of the specimen is preset with a ${60}^0$ angle V-shaped groove with a depth of 1mm, whereas the bottom is preset with a ${20}^0$ angle weakening groove with a depth of 1.5mm. The experimental setup is shown in Figure 5a. The distance between the loading hammer and the weakening groove is $L=24.5$mm, and the hammer loads slowly and vertically downwards, which can be considered as quasi-static loading. The focus of this simulation is on the crack mode. The cracking position of the structure is far away from the load application region, and its deformation is only related to the resultant force and moment of the load. According to the Saint--Venant principle, the specific distribution of the load only affects the stress distribution around the load application region. The error caused by this approximation is confined near the loading boundary. Therefore, considering the requirements of computational efficiency, we scaled the length dimension according to the equivalent of the resultant moment at the groove position. On the other hand, because the ratio of the width to the thickness of the experimental specimen is 10:1, we adopted the plane strain assumption to simplify the three-dimensional structure to a two-dimensional numerical specimen. The dimensions of the numerical specimen and the boundary conditions are shown in Figure 5b.

A grid of the geometric model is shown in Figure 6a. The model was discretized using four-node quadrilateral elements with a total of 2093 nodes and 2009 elements. The grid size was 0.1mm in the region delimited by $y\in [0,4],x\in [1.2,4]$ and uniformly transitioned outside of this region from 0.1mm along the x direction to 0.4mm at the two boundaries. Figure 6b shows the value distribution of the Morphing function $w\left(\mathit{x}\right)$. From Equations (26) and (27), it can be seen that the model partitions are determined by the values of $w\left(\mathit{x}\right)$ and $w\left(\mathit{x}\right)v\left(\mathit{x}\right)$. We preset $v\left(\mathit{x}\right)=1$ for the whole body $\mathsf{\Omega }$, that is, the model partition was only affected by the value distribution of $w\left(\mathit{x}\right)$. Thus, in Figure 6b, the red region corresponds to the pure PD model ($w\left(\mathit{x}\right)=1$), the blue region corresponds to the pure CCM model ($w\left(\mathit{x}\right)=0$), and the transition region is the coupled model ($0<w\left(\mathit{x}\right)<1$).

Because the bond-based PD model was used, Poisson’s ratio for plane strain problems was fixed at 1/4 [4,26]. The horizon $\delta$ of the PD model was set to 0.3mm and the displacement increment was $\Delta \overline{\mathit{u}}=0.018\mathrm{mm}$. At 800°C, the Young’s modulus of the GH4099 superalloy is 147GPa, and its thermal expansion coefficient is $15.1\times {10}^{-6}/°\mathrm{C}$. The yield ($s_Y$) and critical ($s_C$) stretch were set to 0.04 and 0.1 respectively. Figure 7a shows the damage contour at the 15th loading step. Owing to the ductility of the superalloy material at 800°C, the initial crack propagation direction is affected by the V-groove at the top of the structure and has a certain angle with the vertical direction, which is different from that of a brittle fracture. As the loading progresses, the crack continues to propagate along the initial direction. By the 20th loading step (see Figure 7b), the crack starts to deflect toward the bottom weakening groove. By the 40th loading step (see Figure 7c), the crack has propagated to the vicinity of the tip of the weakening groove, then the crack propagation direction becomes parallel to the hypotenuse of the weakening groove. In this state, the structure is difficult to completely fracture. Thus, after the 50th loading step (see Figure 7d), the crack does not continue to grow as loading increases.

At 900°C, the Young’s modulus of the GH4099 superalloy is 121GPa and its thermal expansion coefficient is $15.3\times {10}^{-6}/°\mathrm{C}$. The yield ($s_Y$) and critical ($s_C$) stretch were set to 0.04 and 0.12, respectively. The simulation results obtained under a 900°C environment are shown in Figure 8. The whole fracture process is similar to that at 800°C. The direction of crack propagation changes from the 16th step to the 24th step, as shown in Figure 8a,b. Afterwards, the crack propagates to the vicinity of the tip of the weakening groove and the path becomes parallel to the hypotenuse of the weakening groove by the 48th step. After the 54th loading step, the crack stops growing. However, the crack propagation direction deflects in a more pronounced way compared with the results obtained at 800°C. Figure 9 shows a photomicrograph of the GH4099 superalloy structure at 900°C. An arc-shaped crack path can be seen in the figure, which is consistent with the simulated crack path.

To further verify the effectiveness of the proposed model for simulating ductile fracture at high temperatures, the force-displacement curves obtained from the simulations and experiments at 800°C and 900°C were compared, and the results are shown in Figure 10. Unlike in brittle fractures, the force does not quickly decrease to zero after reaching the peak value. The simulation results were in good agreement with the experimental results. In both the simulations and experiments, it can be seen that, at the end of the curve, the force does not decrease with an increase of displacement but remains constant. This indicates that in the final stage of the failure process of superalloy structures, the structure does not further fracture as the loading progresses. This coincides with the simulated final crack mode. By comparing the force-displacement curves of 800°C and 900°C, we can see that the peak force at 900°C is significantly lower than that at 800°Cowing to the softening effect of the high-temperature environment on the alloy. Additionally, the displacements of both crack initiation and final break decrease, which shows that superalloy structures become more prone to failure as temperature rises.

5.2. Example 2: Four-point bending beam

This example aims to illustrate the effect of the temperature field and plastic deformation on crack propagation paths. The geometry and boundary conditions for this numerical example are shown in Figure 11a. The grid of this example is shown in Figure 11b, and the grid size in the PD region was set to 2mm. Figure 11c shows the value of the morphing function $w\left(\mathit{x}\right)$, which indicates the scope of the different models. The fracture process of the four-point bending beam was simulated under two different temperature conditions, namely a uniform temperature field of 500°C (Figure 12a) and a non-uniform temperature field transitioning from 0°C to 1000°C (Figure 12b). This setting ensures that the specimen is always at 500°C in the vertical direction of the notch of the beam.

To consider the influence of plastic deformation and temperature on the crack path of the four-point bending beam, four numerical tests with the same Young’s modulus and thermal expansion coefficient were conducted. The parameter settings for the four numerical tests are listed in Table 1. For all numerical tests, the horizon $\delta$ was set to three times the grid size of the PD region. The example assumes a plane-stress condition, so Poisson’s ratio is fixed at 1/3.

First, we compared the crack paths of the four-point bending beam under different temperature conditions. Figure 13 shows the final crack paths of cases 1 and 3. These two cases had the same parameter settings except for temperature, but nonetheless exhibited different crack modes. It can be seen that under the uniform temperature field (see Figure 13a), the crack propagates along the vertical direction of the notch of the beam, whereas it propagates along a different path under the non-uniform temperature field (see Figure 13b). This is because there is a transverse temperature gradient in the non-uniform temperature field, which leads to non-uniform thermal expansion at the left and right ends of the beam, resulting in the deflection of the crack.

In cases 2, 3, and 4, the numerical simulations are performed under the non-uniform temperature field, but plastic deformation was only allowed in cases 3 and 4. Thus, in case 2, a bond breaks immediately after its deformation reaches the elastic limit (which is considered a brittle fracture). In cases 3 and 4, an admissible plastic elongation of $\overline{s}=0.005$ and $\overline{s}=0.01$ is introduced after the deformation of a single bond reaches the elastic limit, respectively. The elastic limit of the bond was the same for these three cases. Figure 14 shows the final crack paths of the four-point bending beam for cases 2, 3 and 4. By comparing the crack paths of these three cases, it can be seen that the angle of the crack path in case 2 deviates the most from the vertical direction under the non-uniform temperature field(see Figure 14a). In addition, the deflection angle of the crack path decreases with the increase of the admissible plastic deformation of the bond.

We also analyzed the plastic deformation zone of the crack tip during crack propagation in cases 3 and 4. In Figure 15, subfigures (a), (b) and (c) show the plastic deformation zones around the crack tip at the 10th, 15th and 25th loading steps, respectively, corresponding to region I in Figure 14. On the other hand, subfigures (d), (e), and (f) show the plastic deformation zones around the crack tip at the 10th, 15th and 25th loading steps, corresponding to region II in Figure 14. The color indicates the value of the non-local plastic deformation calculated by Equation (6). It can be seen that, compared with case 3, the crack tip in case 4 is accompanied by a greater nonlocal plastic deformation during the crack propagation. This greater nonlocal plastic deformation around the crack tip weakens the influence of the temperature gradient on the crack path. Therefore, in case 4, the angle at which the crack path deviates from the vertical direction is smaller than that in case 3.

6. Conclusion

We proposed a plastic-fracture PD model for fracture simulation in ductile materials that can accurately simulate plastic deformation and the fracture process under different temperature conditions. We simulated the plastic-fracture process of a GH4099 superalloy structure at 800°C and 900°C. The simulation results show that the proposed model predicts an arc-shaped crack path, which was observed in experiments. The simulated fracture modes and force–displacement curves were in good agreement with our experimental results, demonstrating the effectiveness of the proposed model for simulating ductile fracture at high temperatures. We also carried out numerical tests of a four-point bending beam under uniform and non-uniform temperature fields. The effects of temperature and plastic deformation on the crack path are discussed. The results show that the plastic deformation zone around the crack tip weakens the impact of temperature on the crack deflection.