Stable crack propagation in steel at 1173K: Experimental investigation and simulation using 3D cohesive elements in large-displacements

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers. L’archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d’enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.


Introduction
In the event of a serious accident on a pressurized water reactor (PWR) involving fusion of the fuel core and internal structural elements, the bottom head is subjected to significant thermal and mechanical solicitations. The TMI-II accident of March 28th, 1979 at Three Mile Island, USA is a reference in terms of serious accident scenarios due to the data collected during the event. During this accident, 20 tons of highly calorific corium, combining molten metal and oxide, were relocated in the bottom head. This corium bath led to a 1300 K temperature increase of the bottom head in 30 min, with a 1 m-diameter hot spot. During the accident, the pressure on the reactor vessel was 10 MPa. However, in spite of numerical predictions to the contrary, no cracking occurred in the vessel.
In the hypothetical case where the reactor vessel would suffer cracking, it is necessary to predict the time of occurrence of a crack, its size and its location in order to determine an out-of-vessel accident management strategy within the confinement system. Experimental programs [1][2][3] have been carried out in order to characterize these parameters.
In this paper, our objective is to study more particularly the initiation and quasi-static (or stable) Mode I crack propagation in pressure vessel steel (16MND5) at 1173 K [4]. In the course of this work, cracking tests were carried out on sidegrooved CT specimens. We propose a numerical simulation of these tests using a 3D finite element model and the ABAQUS 6.6-1 analysis code [5]. The crack's initiation and propagation are modeled using cohesive surface elements under large-displacement and large-strain assumptions.
The cohesive zone concept was proposed independently by Barenblatt [6] and Dugdale [7] in order to regularize the crack-tip singularity predicted by linear elastic fracture mechanics [8]. Barenblatt reckoned that near the crack's tip the interatomic forces which tend to bring the two surfaces being pulled apart back together cannot be ignored in the analysis. Dugdale introduced an expanded crack-tip zone in order to limit the development of stresses to a few times the yield stress.
Thus, the role of cohesive zone models is to reproduce the mechanical effects of the fracture process zone at the crack's tip. The zone is modeled using two surfaces which, initially, coincide. During the loading, the two surfaces come apart until complete decohesion. Then, they simulate the lips of the crack. A traction-separation law governs the traction between the two surfaces as a function of their separation and, thus, models the fracture phenomena associated with the material being studied.
The application of the cohesive model (CM) to the ductile fracture of metals was initiated primarily by Tvergaard and Hutchinson [21], Siegmund and Brocks [22] and Tvergaard [23]. The traction-separation laws they proposed are appropriate to the phenomena of nucleation, growth and coalescence of voids [24] which are encountered in ductile fracture.
Only few studies have been carried out using cohesive elements in a 3D formulation [25,26] and under large-strain assumptions [27][28][29]. Chen et al. [27], however, showed in an analysis of the tearing of CT specimens that a plane strain simulation is unable to reproduce the triaxiality state at the crack's tip which can be obtained near the middle section of the 3D model using an identical traction-separation law.
Consequently, the two formulations lead to different crack propagation predictions. Therefore, a 3D analysis is necessary. Besides, the high ductility of our material at high temperatures requires the use of a large-displacement, large-strain formulation.
It is commonly admitted that the energy and the cohesive stress are the only relevant parameters which characterize crack propagation in a ductile material. A few authors have cast doubt on this assumption [30,31]. We will show that in our case of a viscoplastic, highly ductile material at 1173 K the form of the traction-separation law affects the propagation significantly.
Section 2 presents the tests on CT specimens and describes the material, the geometry of the specimens, the testing apparatus and the results. Section 3 introduces the cohesive zone model. The large-displacement 3D formulation of the cohesive elements is presented in detail. We also discuss the modification of the form of Tvergaard and Hutchinson's traction-separation law with respect to the viscous behavior of the material.  The modeling of the tests is described in Section 4. Section 5 presents the optimized fitting results of the traction-separation law in order to reproduce the tests. Based on these results, different assumptions concerning the simulation options and the traction-separation law (the need to introduce some viscosity, the dependence on the degree of triaxiality) are validated.
Section 6 proves that in the case of tearing tests of a CT20 specimen at high temperature ð1173 KÞ the form of the traction-separation law has significant impact, not only locally in the fracture process zone, but also on the global response of the tearing tests.

The material and test setup
The material being studied is a ''Krakatoa"-grade 16MND5 steel whose chemical analysis is presented in Table 1. The specimen had an unconventional geometry based on a CT20 (see Fig. 1): some dimensions were increased in order to reduce necking at high temperature. Lateral grooves were made in order to approximate a plane strain state. A precrack whose profile is described in Table 2 was obtained by fatigue.
The tests were carried out under primary vacuum in a furnace consisting of a water-cooled dual-wall chamber. The temperature was controlled through a B-type thermocouple and its homogeneity measured with three K-type thermocouples. The temperature loading was first applied at a rate of 1200 K h À1 up to 1173 K. Throughout that phase, a constant force was maintained to compensate for depressurization. The actual tearing test was performed under prescribed displacement at a rate of 5 mm min À1 in order to minimize load relaxation. The displacement set point was ended at 15 mm of the actuator's range. The load and displacement were measured using the sensors of the testing device. The crack's propagation was determined using an ANS crack monitoring device. This system works by injecting a current into the specimen through welded wires in order to measure the potential difference between the two sides of the defect. A crack's propagation leads to a decrease in section characterized by an increased resistance and, therefore, an increased potential difference [32]. This method was calibrated during preliminary tests performed at different final extensions of the actuator. Besides, the signal thus obtained is noisy and must be post-processed using linear smoothing. In our case, this smoothing was achieved over a 2 mm crack increment. Consequently, the results at both the initiation and the beginning of the crack's path were extrapolated.
The tests were performed twice. The two results were similar in terms of load-displacement and crack's progression curves. Fig. 2 illustrates the results of one of the tests. Fig. 3 shows the deformed shape and the crack's profile at the end of the test.

The cohesive zone model
The cohesive element and its traction-separation law were implemented as a User Element Fortran subroutine into the ABAQUS/Standard analysis code.

Principle of virtual work under large-displacements
A solid subjected at time t to surface loads F over @B 1 text and prescribed displacements u imp over @B 2 text is traversed by a cohesive surface S tcoh (Fig. 4). Neglecting the terms relative to volume forces and inertia, and using the configuration of the solid at time t as the reference, one can express the principle of virtual work as: where d denotes the separation in the cohesive zone, D the strain rate tensor and T the Cauchy stress tensor.

Formulation of the surface cohesive elements under large-displacements and large strains
The isoparametric bilinear interface element, initially with zero thickness, has eight nodes and four integration points. The nodes are divided into two 4-node surfaces representing a portion of the lower and upper lips of the crack undergoing decohesion.
The movement of the upper and lower parts of the bulk induces a displacement of the nodes of the interface element. Only the relative displacements which represent the three fracture modes are taken into account in the formulation. Here, we are focusing on Mode I alone (opening mode), i.e. the mode which is solicited in the tests. These displacements generate a traction force which prevents their development. This traction is governed by the traction-separation law described in the following section.
In a large-displacement, large-strain formulation, the choice of the reference surface is paramount in characterizing the separation mode of the element [33]. In order to define the local coordinate system of the separation, Ortiz and Pandolfi [29] proposed to use the midsurface S t calculated from the nodal positions of the upper and lower surfaces. Besides, the largestrain formulation enables one to take into account the area reduction due to necking of the bulk during the crack's progression. Fig. 5 shows the geometry of the element. Nodes 1-4 belong to the lower side of the element while nodes 5-8 belong to the upper side. Initially, the two surfaces coincide. Each node has three translational degrees of freedom.
The calculation of the separation at any point of the element is associated with a local coordinate system at the center of the element's midplane in its deformed configuration (Fig. 5). This reference coordinate system ðn; g; fÞ is deduced from the global system by a rotation R T 3Â3 . f is the normal to this midplane, and n and g are the two tangent vectors: Now we can express the value of the relative separation d 12Â1 of the two surfaces at each of the four Gauss points G, of coordinates ðn g ; g g Þ in the element's local coordinate system: ; ; ; ; Functions N i correspond to the four usual bilinear polynomial interpolation functions expressed in the local coordinate system of the element.
Then, the traction vector t 12Â1 at the Gauss point can be calculated using the traction-separation law C presented in the next section: After linearization of the principle of virtual work (Eq. (1)) and finite element discretization, the calculation of the internal forces is expressed on the n þ 1 configuration at the end of the step. Let J denote the Jacobian which transforms the element's local coordinates (Domain h) into global Cartesian coordinates.
Using a Gauss method for numerical integration, one gets the nodal forces: where S represents the element's area, and x g the integration weight of Gauss point g.
The tangent stiffness matrix K T is: ðn g ; g g ÞB nþ1 ðn g ; g g Þ: ð9Þ In this expression, D is the tangent constitutive relation which links the stress to the strain. The expression of K T given here is a simplified version valid only for small strains and small rotations.
This definition is sufficient for the case being considered in this work, where the cohesive zone remains in the plane of symmetry and does not rotate.

The traction-separation law
The traction-separation law represents, on the macroscale, the dissipative phenomena leading to the complete decohesion of the crack's lips in the fracture process zone at the crack's tip.
In the case of a metal at low temperature, the traction-separation law represents the local processes of void nucleation, growth and coalescence. The form of the traction-separation law is unknown a priori.
Tvergaard and Hutchinson [20,21] proposed the law described in Fig. 6. This law is relatively simple, with four parameters, and consistent with an irreversible deformation of the process zone and perfect plasticity of the surrounding bulk.
The same authors showed that in the case of a nonviscous elastic-plastic material the global form of their law has little influence on the crack's propagation. Only two parameters (the maximum traction t max and the rate of energy dissipated during the decohesion of the two surfaces C 0 ) have an influence.
A good candidate for studying this influence, called the material's ductile toughness (R), was proposed by Turner and Kolednik [34]. By extension of the concept of energy recovery rate [35], they proposed the following expression for the toughness of a ductile material: where U denotes the internal energy, E e the elastic energy, B the thickness of the specimen, E pl the inelastic energy and a the crack's progression. R represents the dissipated energy required for the extension of the cracked surface by an increment Bda. This energy is divided into two contributions: an energy dE pl da dissipated inelastically in the bulk and a surface energy C 0 dissipated during the creation of the two lips of the crack. The latter contribution corresponds to the parameter C 0 of the traction-separation law. While these two contributions can be easily separated by modeling a crack with cohesive elements, their experimental characterization has scarcely been investigated [36,37].
Let us consider a constitutive relation of the elastic-plastic bulk with any type of strain hardening. For a constant C 0 , the larger the maximum traction t max relative to the yield strength of the bulk, the larger the parameter R because the rate of energy dissipated in the bulk through plasticity increases [30]. For a constant t max , the larger the surface energy C 0 , the larger R, both intrinsically and because the size of the fracture process zone increases and causes an increase in the plastic energy dissipated.
In our investigation, the temperature is very high ð1173 KÞ. Viscous effects in the bulk material become important. Hence, the localization process leading to fracture is time-dependant.
If the loading rate is slow, the damage process leading to fracture includes diffusive mechanisms (grain boundary cavitation and grain boundary diffusion). The traction-separation law is used in this paper to simplify the localization process which leads to fracture. Consequently, the traction-separation law should be described by a viscous model.
If it is fast, as in this experiment (the experimental loading pin velocity is set constant at 5 mm min À1 ), the fracture is of plastic ductile type depending on the strain rate effect [38]. Hence, the traction-separation law should be simply rate dependant elastoplastic curves. However, the only available tests are constant strain rate (loading pin velocity is constant). Therefore, the cohesive zone model will be elastic-plastic and independent of the global opening rate, which is supposed to be constant during the whole test.
Near the crack tip, in the first stage of the propagation, the strain rates increase and results in the strengthening of the properties of the bulk due to the singularity. However, this strengthening is attenuated by the damage in the zone surrounding the crack's tip. In order to account for this phenomenon simply, Tvergaard and Hutchinson's constitutive law was modified by adding a hardening slope, as shown in Fig. 6b.
Besides, the cracking mode of the CT specimen is Mode I. Therefore, mixed modes were not implemented. Spurious displacements in mixed mode are controlled through a simple elastic stiffness and do not affect the response in Mode I.
The interface is considered to be ruptured once all cohesion has been lost. Then, the separation is greater than d c .

Mesh and boundary conditions
The mesh of the CT specimen, which consists of about 46,000 elements, is shown in Fig. 7. The elements which constitute the bulk are brick elements with reduced integration (eight nodes and one integration point). The mesh size transition near the specimen's fillets is provided through wedge elements (six nodes and two integration points).
The elements of the cohesive zone are laid out appropriately in order to follow the profile of the precrack. The surface is meshed regularly with 43 elements in the direction of the propagation by 20 elements in the direction of the front (860 elements). Two opposite nodes near the central position of the precrack are fixed in the direction of the propagation ðU 1 Þ.
The pins are modeled with rigid shell elements with their reference points located along the axis of the pin. Rotations of the pins other than about their axis ðU 2 Þ and translation about this axis are not allowed. The displacement uðtÞ of the axis of each pin in the traction direction ðU 3 Þ is set at 2.5 mm min À1 . Thus, the midplane of the cohesive zone remains orthogonal to ðU 3 Þ and undergoes no rotation (see Section 3.2).

The calculation procedure
A quasi-static calculation was performed, but the loading rate was respected in order to account properly for the viscous phenomena associated with the material.
The subintegration associated with the elements of the bulk resulted in hourglass displacement modes [39] around the elements adjacent to the cohesive zone. These displacement modes can perturb the propagation of the crack because of the formation of grooves whose amplitude increases with the loading. We chose to stabilize these modes using additional stiffnesses.

The constitutive relation of the bulk
Because of the temperature ð1173 KÞ, the material's behavior is highly dependent on the strain rate. This dependence was taken into account by introducing, at different logarithmic strain rates, a scattered pattern of elastic-plastic curves with strain hardening. The code used linear interpolation to define the behavior at intermediate strain rates.
The Young's modulus E of the material is independent of the strain rate and equal to 23,466 MPa. The elastic-plastic characteristics used in the analysis are given in Table 3 and represented graphically in Fig. 8.

Constitutive law of the interface
For the purpose of the simulations, certain parameters characterizing the form of the traction-separation law, given in Fig. 6b, were assigned.
Parameter d 1 was chosen such that: That was because some authors showed the limits of validity of a traction-separation law with zero initial traction [31,40]. The cracking behavior predicted by cohesive zone models with elastic initial response is mesh-dependent and, thus, different from that predicted by models with infinite initial stiffness [31]. We cannot avoid using this type of law because of the solver used by ABAQUS (Newton Raphson). However, we minimized its influence by choosing the highest possible stiffness with respect to the stiffness of the bulk (convergence limit). Besides, forcing the crack to travel through a single cohesive zone reduces the undesirable effects of this type of law.
The value of d 2 was chosen as proposed by Cornec et al. [25]: Under these conditions, the value of d c can be deduced from the choice of the three other parameters t 1 ; t 2 ; C 0 : Now let us study the influence of parameters t 1 and t 2 on the simulation of the rupture.

Simulation of the tests
This section presents the best simulation of the tests, which corresponds to the following choice of the parameters: C 0 ¼ 47; 000 J=m 2 ; t 1 ¼ 164:52 MPa; t 2 ¼ 197:43 MPa. The corresponding traction-separation law is shown in Fig. 9.

Comparison with the test results
The global load-displacement and crack propagation-displacement curves are in very good agreement with the experimental results (Fig. 10). The difference between the simulation and the test at the beginning of the crack's propagation is due to the simplification of the experimental results (see Section 2). in the creation of the crack's lips in the cohesive zone represents only 9.3% of the internal energy at the end of the calculation. Besides, one can observe that the energy used to control the hourglass modes is negligible compared to the energy dissipated in the cohesive zone (%0.4% of the internal energy at the end of the calculation), which justifies the use of this mode of control in the case of our simulation. Fig. 12 shows the distribution of the energy dissipation rate R during the propagation. The numerical calculation was carried out by fitting the curves of the inelastic energy dissipated in volume U UN and the energy dissipated in the cohesive zone U CZM as functions of the average progression of the crack a. The width of the crack's front used was approximated by the initial thickness of the specimen B i (25 mm in the useful zone):

Results in terms of energy: validation of the control of the displacement hourglass modes
One can verify that the rate of energy dissipated in the cohesive zone corresponds to C 0 in the stationary part of the propagation.  The separation velocity of the crack's lips was calculated as the average over the Gauss points located along a line in the propagation direction (U 1 Þ, near the median position of the specimen's width. Only Gauss points including a separation between dc 4 < d < d c were taken into account in the calculation. Fig. 13 shows that the separation velocity remained constant throughout the test (2 Â 10 À2 < V < 4 Â 10 À2 Þ. Therefore, it was unnecessary to use a viscous traction-separation to model this test.

Dependence of the traction-separation law on the triaxiality rate
The local process of nucleation, growth and coalescence of cavities which is specific to ductile rupture has been shown to be dependent on the stress state [41]. An appropriate way to detect this phenomenon is to measure the triaxiality rate H, which can be calculated from the average stress r h and Von Mises' equivalent stress r eq : In Siegmund and Brocks' works [22], the traction-separation law of the cohesive zone depends on the triaxiality rate, which is assumed to be equal to that of the adjacent element layer in the bulk. Then, the parameters (surface energy dissi-   pation rate and maximum stress) are modified. Based on a 2D plane strain analysis, Siegmund and Brocks showed that in the case of a CT specimen (high-constraint specimen) it is not necessary to take this dependence into account. This result is no longer valid for MT specimens (low-constraint specimens).
However, Chen et al. [27] showed on a 3D analysis of CT10 specimens that the variation of the triaxiality rate along the crack's front throughout the thickness must be taken into account in order for the crack's propagation, and particularly the tunneling effect, to be represented correctly.
We calculated the triaxiality rates for our specimen's geometry (side-grooved CT 25). Fig. 14 shows the maximum triaxiality rate calculated near the Gauss points of the upper layer of bulk elements adjacent to the cohesive zone. Fig. 15 illustrates the position of the corresponding Gauss points and the position of the crack's front in the ðU1; U2Þ plane (Fig. 7). The triaxiality rate is constant over 80% of the width of our specimen, compared to 50% in the case of the CT10 specimen calculated by Chen et al. Therefore, we chose not to make the traction-separation law dependent on the triaxiality rate.

Dependence on the form of the traction-separation law
As already mentioned in Section 3.3, it is usually agreed that for a ductile material only the surface energy recovery rate C 0 and the maximum stress t 2 have an influence on the crack's propagation. However, for the simulation of the tests, it was necessary to add an additional slope to the traction-separation law developed by Tvergaard and Hutchinson. Therefore, we studied the influence of this slope on the propagation results for an identical constitutive relation of the bulk, keeping the values of C 0 and t 2 the same. The results are given for 3 values of t 1 . The corresponding traction-separation laws are shown in Fig. 16.    Table 4 show that the influence of the second slope is anything but negligible. The crack's initiation is not affected, but the propagation velocity is influenced by this parameter. The larger the slope, the higher the propagation velocity. This leads to a decrease in the global stiffness of the specimen. Besides, Fig. 17 illustrates the influence of the slope on the evolution of the toughness (R) during the propagation. For a fixed surface energy dissipation rate, the larger the slope of the traction-separation law, the lesser the energy rate dissipated in the bulk. This results in a decrease in the toughness. Besides, one can observe a decrease in the toughness during the crack's propagation. The smaller the slope of the traction-separation law, the sharper this decrease.

Local influence
The process zone is defined as the zone in which the separation is between d 1 and d c . Fig. 19 shows the evolution of the traction load in the median position of the width of the cohesive zone. This evolution is plotted for a similar crack's progression of about 4 mm. Locally, one can observe a great variation in the distribution of the traction stress. Besides, the size of the fracture process zone also varies. The smaller the traction load t 1 , the larger the process zone. One can see that the size of the fracture process zone is not negligible compared to the size of the crack. This can explain the great influence of the form of the traction-separation law on the crack's propagation and on the load-displacement curve. Fig. 20 illustrates the influence of the slope on the triaxiality rate (calculated as in Section 5.4). One can see that the form of the law has no influence on the maximum value of the triaxiality rate. The triaxiality rate is small at the crack's tip, rises sharply upstream, then remains constant over the remaining length of the fracture process zone. Therefore, the constant part is longer for a small traction load t 1 . Fig. 21 shows the results in terms of average stress and equivalent stress. The average stress is highly dependent on the value of t 1 and follows the same evolution as the traction load calculated near the cohesive zone (Fig. 19). Indeed, it is the maximum principal stress at the bulk level which governs the traction load near the cohesive zone. The equivalent stress is maximum near the crack's front d ¼ d c , which explains that the triaxiality rate is small at that location. The maximum value of the equivalent stress is rather dependent on the form of the traction-separation law. However, in the rest of the fracture process zone, the influence of the slope is negligible.

Conclusion
This paper proposes a simulation of the tearing of a CT specimen at 1173 K. The crack is modeled using cohesive elements. The specificity of such an analysis lies in the viscoplastic nature of the material and the large strains induced by the temperature of the test. These specificities are taken into account in the cohesive elements at different stages.
The formulation of the elements is three-dimensional. We defined a local coordinate system based on the midplane calculated in the deformed configuration of the element [29]. Indeed, the choice of the coordinate system for the decomposition of the separation is very important. It defines the rupture mode of the element. Thus, due to the large strains, the coordinate system must follow the deformed configuration of the element and take into account the physical reality of the fracture phenomenon. Finally, the formulation takes into account the change in geometry induced by the large strains by calculating the residual and the stiffness matrix in the deformed configuration of the element and by taking the variation of its area into account.
The traction-separation law represents, in a macroscopic way, the local dissipative phenomena leading to the creation of the crack's lips. We used a traction-separation law of the Tvergaard and Hutchinson type [20], modified by introducing a hardening slope in place of the plateau. Indeed, the mechanical properties of the material strengthen as one gets closer to the crack's front because the strain rate increases due to the singularity. However, this phenomenon is lessened by the local damage near the singularity. This approach is validated by the fact that the separation velocity of the crack's lips, on average within the process zone, remained constant throughout the test.
The choice of the form of the traction-separation law is often arbitrary. Indeed, in the case of ductile fracture, it is commonly accepted that only the surface energy recovery rate and maximum traction parameters have an influence on the crack's propagation. In the case of our work, that assumption was not verified. We showed that the form of the traction law plays a significant role.
For given values of the surface energy recovery rate and maximum traction load, the stiffening slope influences the crack's propagation velocity. If one measures the toughness in the Turner and Kolednik sense [34], the larger the slope of the traction-separation law, the less tough the specimen. Therefore, the larger the slope, the higher the propagation velocity and the lower the overall stiffness of the specimen. lines : σ eq linespoints : σ h