Numerical modelling strategies using implicit and explicit methods to simulate quasi-static and dynamic three-points bend fracture tests of a ductile steel

Three-point bend fracture tests have been conducted at different loading rates with a quadratic martensitic steel. The failure energy has been found to increase with loading rate. To get insights in this increase a numerical investigation has been undertaken with different strategies using ABAQUS and IMPETUS softwares in order to address quasi-static and dynamic loading conditions. Simulations were conducted with the ABAQUS software in order to carry out a comparative analysis of both implicit and explicit approaches. In addition to standard Finite Element Method (FEM) applied to quasi-static and dynamic conditions, the eXtended-Finite Element Method (X-FEM) was applied to quasistatic conditions. In both approaches, implicit and explicit, crack initiation and propagation were governed by a critical plastic strain threshold combined with a displacement-based damage evolution criterion. Simulations conducted with the IMPETUS software use an explicit approach and second order elements for both quasi-static and dynamic loading conditions. A node-splitting method using an energy-based damage criterion was employed to simulate the crack initiation and propagation. Experimental data and numerical results have been compared, allowing to determine the ability of these two softwares to simulate accurately three-point bend fracture tests.


Introduction
The influence of loading rate on the fracture strength of a quadratic martensitic steel has been investigated with threepoints bend fracture tests. To interpret the experimental data, a numerical investigation has been undertaken with different computational approaches in order to deal with quasi-static and dynamic loading conditions. A special attention was given to the influence of numerical methods on the results, considering that standard constitutive laws and "off-the-shelf" fracture criteria have been retained for the studied ductile martensitic steel. In that perspective, two softwares have been considered: ABAQUS [1] and IMPETUS [2], the former allowing both implicit and explicit approaches not only with standard finite elements but also with X-FEM elements (yet, limited to the implicit approach) [3], the latter, despite a limitation to pure explicit approaches, provides high order fully integrated finite elements.
The objective of this work is to investigate the incidence of implicit and explicit approaches on the computational result in order to implement tailor-made physical models to represent the studied material. One has to keep in mind that such an approach does not pretend to be exhaustive: in particular, the scope of this work is limited to finite element methods and doesn't include more recent mesh-free approaches like Smooth-Particle Hydrodynamics (SPH) for example. As a matter of fact, further additions to this analysis should be considered in a next future.

Three-point bending test Experimental procedure
For both quasi-static and dynamic loading conditions, the 3 points bend specimen has a 60 mm length and a squared section 5 x 5 mm 2 . The 45° V notch has a depth of 2 mm, see Fig. 1a. The span length characterizing the contact was 50 mm. The quasi-static test was realized with an universal testing machine providing the load and displacement history, see Under dynamic loading conditions, a Charpy like experiment was used, see Fig. 2 [4]. The specimen is impacted by a striker bar 10 mm in diameter and 84 mm in length at velocity ranging from 25 to 160 m/s. An insert with a hemispheric shape is placed against the striker to provide smooth loading of the specimen. The specimen is placed against two Hopkinson pressure bars 10 mm in diameter as to record the contact loads. Both striker and Hopkinson pressure bars are made of a swaged 91%W alloy. The striker guided by a Teflon sabot is sent using a gas gun. The impact velocity is first monitored with a two laser beams device placed 40 mm before impact. A specific setup has been developed, in one hand to arrest the Teflon sabot, and in another hand to guide the striker up to impact through continuous sliding of the striker in the Teflon sabot. The velocities of the striker at impact and after impact are deduced from high speed camera observations. The difference between the kinetic energy of the striker before and after impact provides the energy consumed by the round bar Charpy specimen. The kinetic energy of the broken specimen was found to be less than 1% of the total energy absorbed by the specimen. Consequently, the Charpy energy, K, was estimated by: with ms is the mass of the striker, Vi is the velocity of the striker at impact, Vr is the residual velocity of the striker, and S the cross-section area of the three-point bend specimen.

General considerations about computation techniques involved
It is known that quasi-static simulations are generally performed with an implicit approach whereas most dynamic simulations involve an explicit one. The reason is either that implicit calculations applied to dynamic analyzes or explicit calculations applied to quasi-static analyzes leads to unrealistic calculation times.
To differentiate implicit and explicit approaches, it is necessary to start from the principle of virtual work, as expressed in Eqn (2).

(2)
where is the Cauchy stress matrix, b the internal forces resulting vector, T the surface forces vector, ρ the mass density, v the velocity vector and, finally, the acceleration vector.
Once discretized on a mesh, the equation above is solved by means of a numerical method whose general expression is as follows: where Y(t) is the current system state, Y(t+∆t) the system at the later time, F(t,Y) and G(t,Y) Cauchy-continuous functions and α a scalar ranging from zero to one.
On one hand, when α value is zero, the approach is said implicit [5]: the solving of the equation at the next timestep involves both the current state of the system and the later one. This approach, when applied to linear or quasi-linear problems, is not time-dependent, favoring a quick convergence of the calculation. Therefore, despite being harder to implement, the implicit approach is much more efficient to deal with quasi-static analyzes.
On the other hand, when α value is one, the approach is said explicit: the solving of the equation at the next time step involves only the current state of the system, making it easier to implement. However, CFL stability imposes a strong limitation to the time step size as regards to the mesh size and bulk velocity, making this approach quite impractical for most quasi-static applications.
However, in the scope of the quasi-static bending test, if the initial loading prior to breakup is quite linear, the crack initiation and propagation induce strong singularities of the stress-strain fields at the crack tip. Consequently, in such a context, the use of a pure implicit approach is not straightforward as it will lead to convergence problems whereas an explicit modelling strategy doesn't face these limitations.
In the case of the implicit approach, one way to deal with fracture-involved non linearities is to have recourse to eXtended Finite Element Methods (X-FEM) whose purpose is to embed the crack induced strong discontinuity within the element. In that perspective, the nodal degrees of freedom are enriched with additional displacement terms, as expressed in Eqn (4).

(4)
where u represents the displacement at the integration point, NI the shape function and uI the nodal displacement. Functions H(x), Fα(x) and parameters aI and are inherent to the X-FEM method.
Concerning the explicit approach, a dramatic reduction of the calculation time can be obtained by the implementation of a mass-scaling technique which artificially increases the material mass density so that the resulting time step becomes "acceptable". A downside effect of such a method is the generation of inertia effects which can result into a significant amount of an unwanted kinetic energy. Therefore, a special attention must be paid to ensure that this induced kinetic energy stays within certain limits, typically 1% of the total internal energy.

points bending numerical models
The generic computing model used for all quasi-static simulations, regardless of the approach or the software used is illustrated in Fig. 3.

Fig. 3. Generic computing model used for all quasi-static simulations
The nature of the elements used is "approach and software" dependent: -ABAQUS implicit approach is conducted with quadratic fully integrated hexahedral elements, -ABAQUS explicit approach is relying on linear reduced integration hexahedral elements, -IMPETUS explicit approach is conducted with quadratic fully integrated hexahedral elements.
In addition, as illustrated in Fig. 3, a mesh refinement is applied in the notch area.

Material models and crack propagation
In both approaches, implicit and explicit, and with both softwares, ABAQUS and IMPETUS, the material constitutive model is based on a tabulated plasticity flow stress versus plastic strain coupled with a tabulated damage parameter D. Both data sets are extracted from quasi-static uniaxial tests performed on 8 mm diameter and 40 mm length round samples. The damage parameter D induces a progressive loss of material properties, as illustrated in Fig. 4. In ABAQUS, as shown in Fig. 4, the damage parameter D is triggered by a critical plastic strain value, the crack opening occurring when D equals 1. In IMPETUS, this damage parameter is initiated once the plastic energy reaches a given ratio of the plastic energy at failure implemented in the Cockcroft-Latham criterion, as expressed in Eqn. (5) [7]. The crack opening then occurrs when this plastic energy at failure is attained.

(5)
where Dd is the ductile Cockcroft-Latham damage parameter, Wc is the plastic energy at failure, σ1 is the maximum principal stress and is the effective plastic strain.
Then, the modelling of the crack propagation is "approach and software" dependent: -In the case of the implicit X-FEM calculations performed with ABAQUS, the crack opening is modelled through a further (and imposed) enrichment of nodal displacements of the cracked element, -In the case of the explicit (coupled with mass-scaling) calculations performed with ABAQUS, the crack opening is represented by the numerical erosion of the cracked element, -In the case of the explicit (coupled with mass-scaling) calculations performed with IMPETUS, a mass conservative node-splitting technique is used.

General considerations about computation techniques involved
Dynamic calculations involve strong linearities for which a pure explicit approach is much more suitable. As regards to the time scales involved (a few milliseconds depending on the solicitation speed), a mass-scaling is no longer needed. As a consequence, all the dynamic simulations conducted on both softwares use a pure explicit approach.

Dynamic 3 points bending numerical models
The generic computing model used for all dynamic simulations, regardless of the approach or the software used is illustrated in Fig. 5. The nature of the elements is "software" dependent: -ABAQUS calculations are based on linear reduced integration hexahedral elements, -IMPETUS calculations are conducted with quadratic fully integrated hexahedral elements.

Material models and crack propagation
With both softwares, ABAQUS and IMPETUS, the material behavior is represented by means of a Johnson-Cook constitutive model [8] which defines flow stress as: where σy is the Von Mises flow stress, εp is the equivalent plastic strain, is the equivalent plastic strain rate, * is the normalized strain rate, T is the normalized temperature, A is the initial yield strength, B is the hardening parameter, n is the hardening exponent, C is the strain rate parameter, and m is the thermal softening exponent.
This flow-stress model allows considering, to a certain extent, dynamic hardening and thermal softening phenomena. Associated parameters have been fitted from dynamic compression tests conducted on diameter 9 mm height 5 mm compression samples. These parameters are as follows: A=810 MPa ; B = 600 MPa ; C = 0.07 ; n = 0.12 ; m = 1.
The treatment of the crack initiation and growth process is identical to the ones used for quasi-static explicit simulations.

Experimental results
The failure energy of the quadratic martensitic steel as a function of impact velocity Vi is reported in Fig. 6. A significant increase of the failure energy is observed with the impact velocity. At the highest impact velocity of 360 m/s, the failure energy is about ten times the conventional Charpy energy associated with an impact velocity of 10 m/s. Such an increase has been previously reported for FCC-BCC tungsten alloys and a FCC austenitic steel [4,9]. For the FCC-BCC tungsten alloy, this increase has been interpreted from post mortem analyses revealing an initiation process involving multiple cracks [9]. To get insights of the increase of the failure energy, optical and microscopic observations of the fracture specimens have been conducted, see Fig. 7. The failure is ductile involving a similar void coalescence and growth process. The voids are 2 to 5 µm in size disregarding the impact velocity, see Fig. a2, b2 and c2. Such observations reveal that, with the hypothesis that the initiation process involves a single crack, the failure energy will have to be constant.
Based on observations conducted with a FCC-BCC tungsten alloy, it is possible that the increase of the failure energy with the increase of the impact velocity is associated with an initiation failure process involving multiple cracks [8]. The optical fractographies go along with this hypothesis when considering the macroscale failure process. Under quasi-static conditions, the fracture facies is flat when comparing to the machined notch. Under high impact velocity a step is observed between the machined surface and the fracture facies revealing large plastic deformation. Such high level of plastic deformation favor shear lips which is a signature of high failure energy measured at high impact velocity.

Quasi-static simulations
Experimental and numerical results obtained with both softwares are confronted to the experimental response in Fig. 8. Explicit calculations show an acceptable accordance with experience. The fitted parameters for the damage models are as follows: -For the ABAQUS calculation, the critical plastic equivalent strain is 21%, -For the IMPETUS calculation, the Cockcroft Latham criterion is fitted to 750.10 6 J/m 3 A direct comparison of these criteria and, hence, of the two softwares is not trivial but it appears that the critical plastic strain at which damage is initiated is quite equivalent in both softwares: this critical strain εc is found to be around 18% in IMPETUS, so quite close from the 21% used in ABAQUS. In addition, further analyzes are currently in progress in order to assess and compare local plastic strains and energies when fracture occurs in the vicinity of the notch in both calculations.
It must be also mentioned that the mass-scaling induced a low level of inertia effects, with a kinetic energy to internal energy ratio remaining under 0.6%. As such, the quasi-static conditions hypothesis is verified.
The implicit X-FEM approach, however, doesn't fit completely the experimental curve: indeed, if damage initiation seems to be triggered quite accurately (as proven by a pretty good fit with the experimental curve at first displacement increments), then the simulation diverges and fails to represent the final breakup. However, it must be kept in mind that this approach has been conducted with off-the-shelf models recognized to smooth stress-strain fields and, therefore, too dissipative in the scope of ductile failure propagation. A potential way of improvement could be the use of smaller step times in order to get a more accurate numerical response with, most likely, a dramatic increase of the computing time. Several authors ( [5], [6]) have developed dedicated methodologies to deal with ductile fracture by means of X-FEM, as regards to the main benefit of this method: the crack profile is no longer mesh-dependent.

Dynamic simulations
Dynamic results are presented in the Fig. 9 which exhibit the fracture energies KV obtained at three different striker impact speeds. The increase of failure energy with impact velocity, as previously observed, is well described by the ABAQUS and IMPETUS fitted models. As expected, associated damage parameters increase with the impact velocity, see Table 1.

Table 1. Fitted ABAQUS and IMPETUS damage parameters for studied QS and dynamic conditions
As aforementioned, a direct comparison between these two criteria is not straightforward and further analyzes are currently in progress, in the scope of these dynamic conditions, to get a more detailed feedback. In particular, local plastic strains and energies will be assessed in the vicinity of the notch for both calculations. However, critical plastic strains at damage initiation have already been compared in both codes for the 82 m/s impact configuration: this critical strain in IMPETUS is found to be 41%, i.e. very close to the implemented ABAQUS value.
Once more, the selected damage models are pure phenomenological ones and can't pretend to represent physical phenomena involved in the observed increase of fracture energy with solicitation speed. In particular, a special attention will have to be paid not only to the strain-rate dependence but also to the stress triaxiality in further investigations.

Conclusion
The failure energy of three-point bend fracture specimens made of a quadratic martensitic steel was found to increase with the increase of the applied displacement rate, from quasi-static loading conditions, 0.002 m/s, to dynamic conditions 82 m/s. This increase was interpreted through an increase of the macroscopic deformation of the notch with displacement rate. To get insights in this failure energy increase, a numerical investigation was undertaken with different strategies using the ABAQUS and IMPETUS softwares in order to address quasi-static and dynamic loading conditions. The former software includes both implicit and explicit approaches whereas the latter is limited to a pure explicit resolution scheme.
Numerical analyzes have been split in two steps, with a first step addressing quasi-static conditions, and a second step dealing with the dynamic loading conditions. For both analyzes, the ductile fracture propagation has been modelled through a damage-driven law, this damage being initiated by a threshold plastic strain value for the ABAQUS approach and by a given plastic energy for the IMPETUS approach.
In the scope of the quasi-static loading condition, an implicit X-FEM approach uses with the ABAQUS software leads to reproduce the damage initiation process. However, the simulation diverges and fails to represent the final breakup due to a smoothing effect. A potential way of improvement could be the use of smaller step times in order to restitute more accurately stress-strain field singularities and, as a consequence, valorize this approach recognized to be not meshdependent. However, the quasi-static loading condition has been correctly reproduced with an explicit approach for both ABAQUS and IMPETUS softwares, with the successful use of a mass-scaling technique to reduce calculation time. Additional numerical simulations are currently in progress to allow for a complete comparison of the two codes in terms of local plastic strains and energy responses.
In the scope of dynamic conditions, both softwares led to satisfying results. More precisely, the increase of fracture energy with loading rate was reproduced with the increase of the plastic-strain for the ABAQUS damage initiation criterion and of the energy-based IMPETUS criterion.
These simulations have demonstrated the ability of these two softwares to model three-point bend fracture tests as a function of the displacement rate. Furthermore, these simulations provide perspectives to refine fracture modelling, with the mesh-independent X-FEM method, embedded in the ABAQUS software, and the node-splitting technique associated with high order elements in the IMPETUS software. The latter is recognized to be more energy-conservative than standard numerical-erosion based approaches. In addition, for dynamic loading conditions, additional numerical investigations are required with advanced meshless Smooth-Particle Hydrodynamics (Gamma-SPH) which is a method available with the IMPETUS software. This method will allow, as the ABAQUS X-FEM approach under quasi-static conditions, to provide a mesh independent strategy.