Thermomechanical Peridynamic Modeling for Ductile Fracture

In this paper, we propose a modeling method based on peridynamics for ductile fracture at high temperatures. We use a thermoelastic coupling model combining peridynamics and classical continuum mechanics to limit peridynamics calculations to the failure region of a given structure, thereby reducing computational costs. Additionally, we develop a plastic constitutive model of peridynamic bonds to capture the process of ductile fracture in the structure. Furthermore, we introduce an iterative algorithm for ductile-fracture calculations. We present several numerical examples illustrating the performance of our approach. More specifically, we simulated the fracture processes of a superalloy structure in 800 ℃ and 900 ℃ environments and compared the results with experimental data. Our comparisons show that the crack modes captured by the proposed model are similar to the experimental observations, verfying the validity of the proposed model.


Introduction
Many ductile materials, such as superalloys, are regularly subjected to high temperatures for long periods of time, resulting in irreversible plastic deformation and even ductile fracture. Although several numerical studies on the failure of ductile materials have been conducted [1][2][3], the methods proposed therein are not suitable for predicting crack initiation and propagation processes. This is because that the constitutive models used and the corresponding numerical simulations were based on classical continuum mechanics (CCM). Continuity assumption is one of the basic assumptions of CCM, which requires that the field variables are continuous functions of space and time. The stress divergence term in the equilibrium equation described by CCM requires the second derivative of displacement. However, in the fracture problem, the field variables such as displacement at the crack tip are discontinuous. Therefore, the displacement function at the crack tip is nondifferentiable when CCM is used to deal with the fracture problem, that is, the singularity problem appears. 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 simulations based on PD are gradually becoming apparent [5,6].
PD theroy assumes that interactions still exist between non-contact material points separated by a finite distance. It 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 better suited for simulating fracture processed and do not require complex crack propagation criteria [9].
However, nonlocal integration substantially increases calculation costs, and traditional traction boundary conditions cannot be applied directly in PD models. An effective way to interaction between any point x and a point x in its neighborhood is described by the pairwise force function f . The balance equation at point x can be written as Equation (1), where b is a prescribed body force density and H δ (x) ⊂ Ω is a circular region with x as its center and δ as its radius (see Figure 1). Thus, δ is referred to as the horizon [4]. The pairwise force function f can be defined as Equation (2) f wheref (x , x) is the partial interaction due to the action of point x over point x and, correspondingly,f (x, x ) is the partial interaction due to the action of point x over point x . By introducing plastic deformation and thermal expansion into the bond deformation calculation, the partial interaction can be written as Equation (3) f where ξ = x − x is the relative position vector herein referred to as "bond", η ξ (x , x) = η(x , x) · e ξ denotes the projection of the bond deformation at point x over the bond, η(x , x) = u(x ) − u(x), and e ξ = ξ /|ξ|. In addition, u is the displacement and e ξ is the unit direction vector of the bond. In Equation (3), a(x, ξ)T(x, ξ) represents the deformation of the bond due to thermal expansion [34], where a(x, ξ) is the coefficient of thermal expansion of the bond in the PD model. We defineT(x, ξ) as the effective temperature of the bond. In other words,T is determined by T(x) and T(x ), which are the temperature variations between the current temperatures and the reference temperatures of points x and x , respectively. Thus, we have Equation (4).
Moreover,η ξ =η · e ξ , whereη is the historical plastic elongation of the bond defined by where s Y is the yield stretch. 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 a certain 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 [20]. Moreover, the fracture behavior of materials can be described by bond failure [35]. It is assumed that a bond breaks irreversibly if its stretch s exceeds the critical value s C (generally, s C > s Y ), where s = |ξ+η|−|ξ| |ξ| . After introducing plastic elongation at the bond level, we can calculate an average measure of non-local plastic deformation at a given point as Equation (6) wheres = |ξ+η|−|ξ| |ξ| denotes the plastic stretch of the bond.

Governing Equations
We divide the whole solution region into three non-overlapping subregions, i.e., The PD model is confined to subregion Ω 2 , where plastic deformation and fracture may occur, and the CCM model is employed in subregions Ω 1 to reduce computational costs. The two models transition through subregion Ω m . As shown in Figure 2, the displacement boundary conditionū is applied to part Γ u of ∂Ω, and the traction loadF is applied to part Γ F of ∂Ω. n is a unit vector representing the outer normal, and the entire solution region Ω is affected by the body force b and the steady-state temperature field T(x). The Morphing method is based on a set of governing equations(Equations (7)-(13)) that are valid for any point in Ω. It defines the material properties of the two models as functions that change with position [33]. The governing equations are as follows.

Stiffness/Thermal Modulus Constraint Equations
The governing equations of the hybrid model were presented in Section 3.1. Functions w(x), v(x) and E(x), β(x) jointly guide the transition between the CCM and PD models. w(x) and v(x) determine the division of the hybrid model; thus, how we define E(x) and β(x) is key to successfully couple both models together. In this subsection, we derive the constraint functions for E(x) and β(x) based on the point-by-point equivalency of the energy density.
Plastic deformation is confined to the pure peridynamic region, that is,η ξ ≡ 0 outside this region. Thus, the constitutive equation defined in Equation (13) can be simplified to an elastic version. Let us consider a subregion Ω 0 = {x ∈ Ω : H δ (x) ⊂ Ω}, such that we still have a complete integral region when the PD model is used near the boundary. We can derive the following equations based on a previous study [33]. The following assertions hold true for any point x ∈ Ω 0 : , the model is restricted to the pure CCM model at point x. Then the elastic energy density at point x can be written as Equation (17) W • If and only if • If and only if E(x) = 0, β(x) = 0 and ∃x ∈ H δ (x) such that 0 < w(x) < 1, 0 < v(x) ≤ 1, the hybrid model must be used at point x. Considering Equation (4), the elastic energy density at point x can be written as Equation (19) W Remark 1. W 0 (T) 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(x) and v(x) and only depend on the displacement and temperature field. This implies that the elastic energy densities calculated at the same point with different expressions are equivalent to each other; thus, the corresponding terms of force and temperature must be equal. Therefore, the terms in Equation (17) are equal to the corresponding ones in Equation (18), then we have Equations (20) and (21) Furthermore, we assumed that both the temperature and elastic deformation in a small neighborhood around point x are uniform; thus, we have Equations (22) and (23) Thus, Equations (20) and (21) lead to Equations (24) and (25) Similarly, the terms of Equation (17) are equal to the corresponding ones in Equation (19). By considering Equations (22)-(25) we get Equations (26) and (27) In summary, we have established a thermoelastic hybrid PD-CCM model in Equations (7)- (16), (26) and (27). The partitioning of the hybrid model is determined by the Morphing functions w(x) and v(x), and the definitions of w(x) and v(x) are independent of each other.

Plastic-Fracture Calculations
In this study, the constraint equations of the hybrid model ensure that PD is applied 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 a flowchart of the implicit iterative algorithm used for the plasticfracture simulation. In this algorithm, the coupling and the CCM zones are discretized using continuous finite elements, whereas the pure PD zone is discretized using discrete finite elements, i.e., nodes are not shared between different elements. Discrete elements allow the cracks to grow freely along their boundaries. In fact, we can also discrete the PD zone by continuous elements, then cracks can pass through the element and manifest as damage. The discrete strategy may affect the crack propagation path, but the crack propagation simulated by these two discrete strategies is consistent when appropriate mesh fineness is used. The two discrete strategies are discussed in the reference [36].
The boundary conditions are gradually imposed in N increments. At each incremental step, the elongation of each bond is calculated based on the displacement field obtained by solving the corresponding linear equations. If any bond yields, the residual plastic force must be calculated to update the force vector on 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, we must update the stiffness matrix by removing the contribution of these bonds. Change in the stiffness matrix will inevitably alter the calculated displacement field; therefore, we must recalculate the residual plastic force and solve the equations until the plastic iteration converges and no bonds break at the current step. Subsequently, the step counter is incremented and calculations are repeated until the last load step is completed.

Numerical Examples
In this section, we validate the proposed model by presenting two numerical examples. The first is a two-dimensional simplified model of an actual experiment involving a GH4099 superalloy structure in 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. The numerical examples are calculated by the home-made program based on C++, which is designed in parallel using OpenMP. All calculations were performed on the 12th Gen Intel(R) Core(TM) i7-12700F processor.
It is important to note that the zone in which PD is implemented in the numerical simulations is empirically determined. In previous studies, Wang, Han and Lubineau proposed adaptive coupling modeling techniques based on broken-bond drive [31] and strength drive [32], to update the PD zone and remesh so as to track the crack propagation path. However, For the numerical examples presented in this study, it is feasible to preset a fixed PD zone that can cover the crack path. Therefore, the adaptive coupling modeling and remesh techniques are not considered in the plastic fracture simulations in this study.

Example 1: GH4099 Superalloy Structure
The geometry of the experimental GH4099 superalloy specimen is presented in Figure 4. As shown in Figure 4b, the top of the specimen is preset with a 60 • V-shaped groove with a depth of 1 mm, whereas the bottom was preset with a 20 • angle weakening groove with a depth of 1.5 mm. The experimental setup is shown in Figure 5a. The distance between the loading hammer and the weakening groove is L = 24.5 mm. First, the structures are heated to the target temperatures (800 • C and 900°C respectively) at the rate of 30°C per minute in a high-temperature chamber, and holding temperature for 5 min. Then the hammer as shown in Figure 5a was vertically loaded downward at the speed of 0.5 mm per minute. The measurement tolerance of the high-temperature chamber is ±2°C. And the measurement resolution of hammer loading force is 0.001N. The focus of this simulation was 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 affects only the stress distribution around the load application region. The error caused by this approximation is confined near the loading boundary. Therefore, considering the computational efficiency requirements, we scaled the length dimension according to the equivalent of the resultant moment at the groove position. On the other hand, because the width-to-thickness ratio 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. We selected a small mesh fineness in PD zone to ensure the precision of crack simulations, while a larger mesh fineness was selected in CCM zone to reduce the calculation consumption. The study on mesh convergence of PD-CCM coupling model based on the Morphing method can be found in the reference [16]. The model was discretized using four-node quadrilateral elements into 2093 nodes and 2009 elements. The grid size was 0.1 mm in the region delimited by y ∈ [0, 4], x ∈ [1.2, 4] and uniformly transitioned outside of this region from 0.1 mm along the x direction to 0.4 mm at the two boundaries. Figure 6b shows the value distribution of the Morphing function w(x). From Equations (26) and (27), it can be seen that the model partitions are determined by the values of w(x) and w(x)v(x). We preset v(x) = 1 for the entire body Ω, that is, the model partition was only affected by the value distribution of w(x). Thus, in Figure 6b, the red region corresponds to the pure PD model (w(x) = 1), the blue region corresponds to the pure CCM model (w(x) = 0), and the transition region corresponds to the coupled model (0 < w(x) < 1).
(a) (b) Figure 5.  Because we used a bond-based PD model, Poisson's ratio for plane strain problems was fixed at 1/4 [4,37]. The horizon δ of the PD model was set to 0.3 mm and the displacement increment was ∆ū = 0.018 mm. The material parameters are shown in Table 1, at 800°C, the Young's modulus of the GH4099 superalloy is 147 GPa, and its thermal expansion coefficient is 15.1 × 10 −6°C−1 . 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 at 800°C, the initial crack propagation direction was affected by the V-groove at the top of the structure and had a certain angle with the vertical direction, which was different from that of a brittle fracture. As loading progressed, the crack continued to propagate along the initial direction. By the 20th loading step (see Figure 7b), the crack started to deflect toward the bottom weakening groove. By the 40th loading step (see Figure 7c), the crack had propagated close to the tip of the weakening groove, and the crack propagation direction became parallel to the hypotenuse of the weakening groove. In this state, It is difficult for the structure to fracture completely. Thus, after the 50th loading step (see Figure 7d), the crack did not continue to grow as the loading increased.  At 900°C, the Young's modulus of the GH4099 superalloy is 121 GPa and its thermal expansion coefficient is 15.3 × 10 −6°C−1 . The yield (s Y ) and critical (s C ) stretch were set to 0.04 and 0.12, respectively. The simulation results obtained in a 900°C environment are shown in Figure 8. The entire fracture process is similar to that at 800°C. The direction of crack propagation changed from the 16th step to the 24th step, as shown in Figure 8a,b. Subsequently, the crack propagated close to the tip of the weakening groove, and the path became parallel to the hypotenuse of the weakening groove by the 48th step. After the 54th loading step, the crack stoped growing. However, the crack propagation direction deflected in a more pronounced manner 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 consistent with the simulated crack path can be observed in the figure.
To further verify the effectiveness of the proposed model for simulating ductile fractures at high temperatures, we compared the force-displacement curves obtained from the simulations and experiments at 800°C and 900°C. The results are shown in Figure 10. Unlike in brittle fractures, force did not quickly decrease to zero after reaching its peak value. The simulation results were essentially consentient with the experimental results, and the error of simulated peak forces at 800°C and 900°C relative to experimental results is 0.9 % and 4.9 %, respectively. In both the simulations and experiments, at the end of the curve, the force did not decrease with an increase in displacement but instead remained constant. This indicates that, in the final stage of the failure process of superalloy structures, the structure does not fracture further as loading progresses. This coincides with the simulated final crack mode. By comparing the force-displacement curves obtained at 800°C and 900°C, we can see that the peak force at 900°C was significantly lower than that at 800°C owing to the softening effect of the high-temperature environment on the alloy. Additionally, the displacements of both crack initiation and final break decreased,indicating that superalloy structures become more prone to failure as temperature increases.

Example 2: Four-Point Bending Beam
This example illustrates the effects of the temperature field and plastic deformation on the crack propagation paths. The geometry and boundary conditions for this numerical example are shown in Figure 11a. The grid used in this example is shown in Figure 11b. The grid size in the PD region was set to 2mm. Figure 11c shows the value distribution of the Morphing function w(x), which indicates the scope of the different models. The fracture process of the four-point bending beam was simulated under two different temperature conditions: 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 was always at 500°C in the vertical direction of the beam notch.
To consider the effects of plastic deformation and temperature on the crack path of the four-point bending beam, we conducted four numerical tests with the same Young's modulus and thermal expansion coefficient. The parameter settings for the four numerical tests are listed in Table 2. In all tests, the horizon δ was set to three times the grid size of the PD region. This example assumes a plane-stress condition; therefore, Poisson's ratio was 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 for cases 1 and 3. These two cases had the same parameter settings except for temperature, but nonetheless exhibited different crack patterns. Under the uniform temperature field (see Figure 13a), the crack propagated along the vertical direction of the beam notch. In contrast, it propagated along a different path under the non-uniform temperature field (see Figure 13b). We presented the strain contours before crack initiation of both cases 1 and 3. Figure 14a shows the strain of ε xx under the uniform temperature condition, in which the beam has a symmetrical strain field. It is obvious that the crack will propagate vertically. In the case of non-uniform temperature condition, the transverse temperature gradient leads to the asymmetric strain field (see Figure 14b), and the value of strain ε xx in the right half of the beam is significantly higher than that in the left half. Therefore, the crack will deviate from the center line to the right. In cases 2, 3, and 4, the numerical simulations were performed under the non-uniform temperature field, but plastic deformation was only allowed in cases 3 and 4. Thus, in case 2, bonds broke immediately after their deformation reached the elastic limit (which is considered brittle fracture). In cases 3 and 4, an admissible plastic elongation ofs = 0.005 ands = 0.01 was introduced after the deformation of a single bond reaches the elastic limit, respectively. The elastic limit of the bond was the same in all three cases. Figure 15 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 deviated the most from the vertical direction under the non-uniform temperature field (see Figure 15a). In addition, the deflection angle of the crack path decreased with an increase in 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 16, 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 15. 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 15. Color indicates the value of the non-local plastic deformation calculated using Equation (6). Compared with case 3, the crack tip in case 4 exhibited 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 deviated from the vertical directionwas smaller than that in case 3.  . Nonlocal plastic deformation contours: (a-c) correspond to the 10th, 15th, and 25th loading steps for region I in Figure 15, respectively, (d-f) correspond to the 10th, 15th, and 25th loading steps for region II in Figure 15, respectively. (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this paper).

Conclusions
We propose a PD modeling approach to simulate the plastic deformation and fracture of structures in a thermomechanical coupling field. The proposed model can qualitatively capture the ductile fracture mode of structures. We simulated the plastic-fracture process of a GH4099 superalloy structure at 800°C and 900°C. The simulation results showed that the proposed model predicts an arc-shaped crack path, which is consistent with experimental observations. The simulated crack path and force--displacement curves were essentially consentient with the experimental results, demonstrating the effectiveness of the proposed model for simulating ductile fractures at high temperatures. We also conducted numerical tests on a four-point bending beam under uniform and non-uniform temperature fields, and discussed the effects of temperature and plastic deformation on the crack path. The results show that the plastic deformation zone around the crack tip weakened the influence of temperature on crack deflection.
The proposed model has the advantage of PD in dealing with fracture problems, and the coupling with CCM improves the efficiency and applicability of the model. The PD-CCM coupled model has great potential in the simulation of plastic fracture in thermomechanical coupling field. However, compared with the CCM model, PD constitutive equations still lack experimental verification and verification. Therefore, further research can focus on the experimental verification and verification of PD constitutive equations.