Finite difference modelling for understanding the hydrogen assisted cracking in virtual slow strain rate tensile tests

1. INTRODUCTION

The use of high strength steels is constantly increasing, specially powered by the automotive industry (Necati and Koc, 2014Necati, Ö., Koc, M. (2014). Promises and problems of ultra/advanced high strength steel (U/AHSS) utilization in automotive industry. 7th Automotive Technologies Congress, Bursa, Turkey. ), due to the weight and cost advantages offered by these materials. In the oil and gas industry, high strength steels are applied in pipelines (Mohtadi-Bonab, 2019Mohtadi-Bonab, M.A. (2019). Effects of Different Parameters on Initiation and Propagation of Stress Corrosion Cracks in Pipeline Steels: A Review. Metals 9 (5), 590. https://doi.org/10.3390/met9050590.) and in submerged structural elements for offshore platforms (Krom et al., 1999aKrom, A.H.M., Koers, R.W.J., Bakker, A.D. (1999a). Hydrogen transport near a blunting crack tip. J. Mech. Phys. Solids 47 (4), 971-992. https://doi.org/10.1016/S0022-5096(98)00064-7.; Billingham et al., 2003Billingham, J., Sharp, J.V., Spurrier, J., Kilgallon, P.J. (2003). Review of the performance of high-strength steels used offshore. 1st ed., HSE Books, Sudbury, UK. ). More specifically, they are employed in jack-ups (Sharp et al., 2001Sharp, J.V., Billingham, J., Robinson, M.J. (2001). The risk management of high-strength steels in jack-ups in seawater. Mar. Struct. 14 (4-5), 537-551. https://doi.org/10.1016/S0951-8339(00)00053-8.) and mooring chains (DNVGL-CP-0237, 2018DNVGL-CP-0237 (2018). Offshore Mooring Chain and Accessories. Class Programme. Approval for Manufacturers. DNG GL AS. Available online: https://rules.dnvgl.com/docs/pdf/dnvgl/CP/2018-07/DNVGL-CP-0237.pdf.). Construction steels for mooring offshore structures are classified in “strength grades” according to their tensile properties and named from R3 to R6 (DNVGL-CP-0237, 2018DNVGL-CP-0237 (2018). Offshore Mooring Chain and Accessories. Class Programme. Approval for Manufacturers. DNG GL AS. Available online: https://rules.dnvgl.com/docs/pdf/dnvgl/CP/2018-07/DNVGL-CP-0237.pdf.), where R3 is the lowest resistance grade. Two steels from the same strength grade are considered equivalent for a given mechanical function.

As high strength steels present risk of Hydrogen Embrittlement (HE) when employed in mooring applications, industry standards assess the susceptibility to hydrogen embrittlement by using tests in which a tensile specimen is immersed in synthetic seawater and tensile tested at very low deformation rates (10^-5 s^-1) in the presence of a hydrogenating cathodic potential (DNVGL-OS-E302, 2018DNVGL-OS-E302 (2018). Offshore Mooring Chain, Offshore Standard. DNV GL AS. Available online: https://rules.dnvgl.com/docs/pdf/dnvgl/OS/2018-07/DNVGL-OS-E302.pdf.). These tests are known as Slow Strain Rate Tensile tests (SSRT). Previous work by the authors, comparing a set of environmental conditions ranging from non-hydrogenating to strongly hydrogenating submerged seawater service, showed that SSRT can cause different cracking patterns depending on the studied steel, even for steels classified with the same strength grade (Artola et al., 2018Artola, G., Arredondo, A., Fernández-Calvo, A., Aldazabal, J. (2018). Hydrogen embrittlement susceptibility of R4 and R5 high strength mooring steels in cool and warm seawater. Metals 8 (9), 700. https://doi.org/10.3390/met8090700.).

Table 1 shows an extract from their work in which two mooring chain steels (Sample A and Sample B) were extracted from full scale mooring chain elements belonging to strength grade R4. Both samples, whose microstructure was composed by a quench and tempered matrix with a minor differences between each other in terms of their carbide dispersion (Artola, 2018Artola, G. (2018). Susceptibilidad a la fragilización por hidrógeno de aceros de alta resistencia: comportamiento en ambientes marinos y modelización de patrones de agrietamiento. PhD. Thesis, Tecnun-Universidad de Navarra, Donostia-San Sebastián. http://hdl.handle.net/10171/56315.), showed similar tensile properties both during conventional testing and under strong hydrogen embrittling SSRT conditions (cathodic overprotection of -1200 mV with the steel immersed in synthetic seawater during the test). Nevertheless, the surface cracking pattern they showed in the hydrogenating condition was different, as shown in Fig. 1. Sample A showed a significantly higher number of surface cracks when hydrogen embrittlement was involved in the fracture.

In HE, the hydrogen atoms that are dissolved in the metal interact with its microstructure by means of a set of mechanisms that modify its mechanical behaviour. Thus, the evolution of the distribution of the hydrogen dissolved in the microstructure determines the expected mechanical response for given load and geometry scenario (Toribio, 2011Toribio, J. (2011). Fracture Mechanics approach to stress corrosion cracking of pipeline steels: when hydrogen is the circumstance. In: Integrity of Pipelines Transporting Hydrocarbons. Springer, Dordrecht, Holland, pp. 37-58. https://doi.org/115-126. 10.1007/978-94-007-0588-3_4 ) and, therefore, affects the cracking pattern.

The amount of dissolved hydrogen in steels can be split in two components: a fraction in interstitial solid solution within the iron matrix and the rest trapped in microstructural elements such as second phases, grain boundaries or dislocations, which are known as “hydrogen traps” (Pound, 1998Pound, B.G. (1998). Hydrogen trapping in high strength steels. Acta Mater. 46 (16), 5733-5743. https://doi.org/10.1016/S1359-6454(98)00247-X.; Park et al., 2002Park, D., Maroef, I., Landau, A., Olson, D. (2002). Retained austenite as a hydrogen trap in steel welds. Welding Journal 81(2), 27-35. Available online: http://files.aws.org/wj/supplement/Park2-02.pdf.; Robertson et al., 2015Robertson, I., Sofronis, P., Nagao, A., Martin, M., Wang, S., Gross, D., Nygre, K. (2015). Hydrogen embrittlement understood. Metall. Mater. Trans. B 46, 1085-1103. https://doi.org/10.1007/s11663-015-0325-y.). The former being known as “diffusible hydrogen”, which can be hypothesized to be the main responsible for the HE, and the latter being known as “trapped hydrogen”. A fraction of the hydrogen can shift from “trapped” to “diffusible” and vice versa depending on the boundary conditions and the energy of the trap. Despite it must be kept in mind that there is an open discussion regarding the generalization of the hypotheses (Oudriss et al., 2014Oudriss, A., Fleurentin, A., Courlit, G., Conforto, E., Berziou, C., Rébéré, C., Cohendoz, J.M., Sobrino, J.M., Creus, J., Feaugas, X. (2014). Consequence of the diffusive hydrogen contents on tensile properties of martensitic steel during the desorption at room temperature. Mater. Sci. Eng. A 598, 420-428. https://doi.org/10.1016/j.msea.2014.01.039.; Depover and Verbeken, 2018Depover, T., Verbeken, K. (2018). The detrimental effect of hydrogen at dislocations on the hydrogen embrittlement susceptibility of Fe-C-X alloys: An experimental proof of the HELP mechanism. Int. J. Hydrog. Energy 43 (5), 3050-3061. https://doi.org/10.1016/j.ijhydene.2017.12.109.; Djukic et al., 2019Djukic, M.B., Bakic, G.M., Sijacki Zeravcic, V., Sedmak, A., Rajicic, B. (2019). The synergistic action and interplay of hydrogen embrittlement mechanisms in steels and iron: Localized plasticity and decohesion. Eng. Fract. Mech. 216, 106528. https://doi.org/10.1016/j.engfracmech.2019.106528.), increasing the hydrogen trapping capacity of the steels has been related to HE resistance improvements, in particular cases with the use of carbides (Yamasaki and Bhadeshia, 2006Yamasaki, S., Bhadeshia, H.K.D.H. (2006). M₄C₃ precipitation in Fe-C-Mo-V steels and relationship to hydrogen trapping. Proceedings of the Royal Society A 462 (2072), 43-50. https://doi.org/10.1098/rspa.2006.1688.). Introducing hydrogen diffusion barriers in the steel has also led to HE behaviour improvement (Fielding et al., 2014Fielding, L., Song, E., Han, D., Bhadeshia, H.K.D.H., Suh, D.-W. (2014). Hydrogen diffusion and the percolation of austenite in nanostructured bainitic steel. Proceedings of the Royal Society A 470 (2168), 20140108-1-17. https://doi.org/10.1098/rspa.2014.0108.). Taking this into account, there is a favourable scenario for hypothesizing that the differences in the cracking patterns observed in the SSRT could be related to hydrogen trapping and diffusion behaviour variations between Sample A and Sample B. For this hypothesis to be feasible, a model for HE considering both total and diffusible hydrogen and its transport through the material should be able to reproduce these cracking patterns.

For practical purposes, the diffusive component tends to be modelled using Finite Element Methods (FEM) in a broad range of detail levels, from stress driven diffusion (Oh et al., 2010Oh, C.S., Kim, Y.J., Yoon, K.B. (2010). Coupled analysis of hydrogen transport using ABAQUS. J. Solid. Mech. Mat. Eng. 4 (7), 908-917. https://doi.org/10.1299/jmmp.4.908.), including complex effects such as the presence of traps (Yan et al., 2014Yan, C., Liu, C., Zhang, G. (2014). Simulation of hydrogen diffusion in weld joint of X80 pipeline steel. J. Cent. South Univ. 21, 4432-4437. https://doi.org/10.1007/s11771-014-2445-y.), to coupled mechanical-diffusion models (Di Leo and Anand, 2013Di Leo, C.V., Anand, L. (2013). Hydrogen in metals: A coupled theory for species for species diffusion and large elastic-plastic deformations. Int. J. Plast. 43, 42-69. https://doi.org/10.1016/j.ijplas.2012.11.005.; Barrera et al., 2016Barrera, O., Tarleton, E., Tang, H.W., Cocks, A.C.F. (2016). Modelling the coupling between hydrogen diffusion and the mechanical behavior. Comput. Mater. Sci. 122, 219-228. https://doi.org/10.1016/j.commatsci.2016.05.030.; Díaz et al., 2019Díaz, A., Alegre J.M., Cuesta I., Zhang, Z. (2019). Numerical Study of Hydrogen Influence on Void Growth at Low Triaxialities Considering Transient Effects. Int. J. Mech. Sci. 164, 105176. https://doi.org/10.1016/j.ijmecsci.2019.105176.). Despite their complexity, these FEM approaches are attractive to the HE study and works whose results are applicable to SSRT have been published (Krom et al., 1999bKrom, A.H.M., Maier, H.J., Koers, R.W.J., Bakker, A. (1999b). The effect of strain rate on hydrogen distribution in round tensile specimens. Mater. Sci. and Eng. A. 271 (1-2), 22-30. https://doi.org/10.1016/S0921-5093(99)00276-2.; Miresmaeili et al., 2010Miresmaeili, R., Ogino, M., Nakagawa, T., Kanayama, H. (2010). A coupled elastoplastic- transient hydrogen diffusion analysis to simulate the onset of necking in tension by using the finite element method. Int. J. Hydrog. Energy 35 (3), 1506-1514. https://doi.org/10.1016/j.ijhydene.2009.11.024.).

Despite recent publications have solved transient diffusion-mechanical processes by FEM (Benannoune et al., 2018Benannoune, S., Charles, Y., Mougenot, J., Gaspérini, M. (2018). Numerical simulation of the transient hydrogen trapping process using an analytical approximation of the McNabb and Foster equation. Int. J. Hydrog. Energy 43 (18), 9083-9093. https://doi.org/10.1016/j.ijhydene.2018.03.179.), previous works coupled Finite Difference Methods (FDM) with FEM when modelling non-steady processes (Ohmi et al., 2012Ohmi, T., Yokobori, A.T., Takei, K. (2012). Flow Controllability of Hydrogen Diffusion Driven by Local Stress Field for Steel under Cyclic Loading Condition. Defect and Diffusion Forum 326-328, 626-631. https://doi.org/10.4028/www.scientific.net/DDF.326-328.626.; Sezgin et al., 2019Sezgin, J., Takakuwa, O., Matsunaga, H., Yamabe, J. (2019). Simulation of the effect of internal pressure on the integrity of hydrogen pre-charged BCC and FCC steels in SSRT test conditions. Eng. Fract. Mech. 216, 106505. https://doi.org/10.1016/j.engfracmech.2019.106505.); being the former used for the assessment of hydrogen diffusion and the later for the assessment of the stresses and strains. FDM has been also used by itself to model diffusion in hydrogen embrittlement related studies (Turnbull et al., 1997Turnbull, A., Hutchings, R.B., Ferriss, D.H. (1997). Modelling of thermal desorption of hydrogen from metals. Mater. Sci. Eng. A. 238 (2), 317-328. https://doi.org/10.1016/S0921-5093(97)00426-7.; Park et al., 2016Park, J.H., Oh, M., Kim, S.J. (2016). Effect of bainite in microstructure on hydrogen diffusion and trapping behavior of ferritic steel used for sour service application. J. Mater. Res. 32 (7), 1295-1303. https://doi.org/10.1557/jmr.2016.480.), for modelling stress-strain fields alone (Lin et al., 2019Lin, D., Fang, X., Gan, Y., Zheng, Y. A. (2019). Damped Dynamic Finite Element Difference Approach for Modelling Static Stress-Strain Fields. Pure Appl. Geophys 176, 3851-3865. https://doi.org/10.1007/s00024-019-02207-2.). Analogous methods employing scalar matrices have been used to model structural damage (Palm et al., 2012Palma, J., Silva, D., Andrade, J.M., Almeida, A. (2012). Numerical modeling of hydrogen diffusion in structural steels under cathodic overprotection and its effects on fatigue crack propagation. Mat.-wiss. U. Werkstofftech. 43 (5), 392-398. https://doi.org/10.1002/mawe.201200971.). But no integrated stress-diffusion-strength modelling approach by FDM has been found in the literature.

Other works simulated successfully the diffusion on interstitial atoms in bcc lattices using random walk techniques (Aldazabal and Aldazabal, 2013Aldazabal, J., Aldazabal, I. (2013). Computer Modelling of Atomic Movements/Diffusion on Oxide Dispersion Strengthened (ODS) Materials. EUROMAT2013, Seville, Spain.; San Sebastian et al., 2008San Sebastian, I., Aldazabal, J., Capdevila, C., Garcia-Mateo, C. (2008). Diffusion simulation of Cr-Fe bcc systems at atomic level using a random walk algorithm. Phys. Stat. Solid. A-Appl. Res. 205 (6), 1337-1342. https://doi.org/10.1002/pssa.200778124.) but the simulated time and length scales were extremely small. Fortunately, it is possible to benefit from existing analogies, from the mathematical point of view, between mass diffusion and heat transfer processes. Recent works (Takaishi, 2017Takaishi, T. (2017). Phase field crack growth model with hydrogen embrittlement. Mathematical Analysis of Continuum Mechanics and Industrial Applications. H.I. Springer Nature Singapore Pte. Ltd., pp. 27-34. https://doi.org/10.1007/978-981-10-2633-1_3.; Martínez-Pañeda et al., 2018Martínez-Pañeda, E., Golahmar, A., Niordson, C.F. (2018). A phase field formulation of hydrogen assisted cracking. Comp. Methods Appl. Mech. Eng. 342, 742-761. https://doi.org/10.1016/j.cma.2018.07.021.; Fu and Fang, 2018Fu, L., Fang, H. (2018). Formation criterion of hydrogen-induced cracking in steel based on fracture mechanics. Metals 8 (11), 940. https://doi.org/10.3390/met8110940.; Hüter et al., 2018Hüter, C., Shanthraj, P., McEniry, E., Spatschek, R., Hickel, T., Tehranchi, A., Guo, X., Roters, F. (2018). Multiscale Modelling of Hydrogen Transport and Segregation in Polycrystalline Steels. Metals 8 (6), 430. https://doi.org/10.3390/met8060430.; Anand et al., 2019Anand, L., Mao, Y., Talamini, B. (2019). On modelling fracture of ferritic steels due to hydrogen embrittlement. J. Mech. Phys. Sol. 122, 280-314. https://doi.org/10.1016/j.jmps.2018.09.012.) have successfully applied mathematical formulations such as phase field modelling, molecular dynamics and ab initio calculations to the nucleation and growth of cracks when HE is in place.

The present work proposes an alternative to the approaches above, which bets for simplicity and implements an iterative calculation strategy based on the FDM. The reason behind suggesting a new modelling strategy among the already existing choices, which lead to successful results, is the accessibility. The simulation proposed model can be reproduced and transferred to other applications with very basic programming knowledge and costless resources. In the case that several fields are considered simultaneously, for example stresses and diffusion, it will be necessary to calculate them carefully to simulate its time evolution properly. This model has been applied to simulate SSRT in materials with different hydrogen trapping capacities and to analyse qualitatively if the differences in cracking patterns resemble the previous experimental observations on R4 structural steel used submerged structures.

2. MATERIALS AND METHODS

The finite difference modelling methodology for the SSRT under HE has been divided in three components that have been superimposed to each other at each step of the simulation, taking into account that the main aim of this model is easy implementation:

2.1. Displacement model

 ⌅

The mechanical component of the model has been approached from the perspective of displacements. First, the domain of interest has been discretized in a homogeneous and regular network of MxN points, elements, cells, nodes or pixels. Next, each node has been assigned a scalar called “equivalent displacement” (ED), u_i,j. This scalar represents the total displacement condition of each element, relative to its initial position.

This leads to a U_MxN matrix, whose elements are initially u_i,j = 0. To model the deformation applied in the SSRT, boundary conditions have been imposed on the top (u_0,j ) and bottom (u_M,j ) sides of the domain. A fixed positive value of the ED, corresponding to the macroscopic vertical displacement to be simulated, has been imposed on the upper face, while on the lower side, the value has been set constant and equal to u_M,j = 0. This resembles a test in which the specimen is fixed at the bottom and pulled from the top. To simulate a tensile test, the value of the imposed displacement, ED, was increased progressively as the stresses were calculated. Depending on the incremental values of ED, the test can be reproduced in a coarser or finer way. High deformation rates will assume that the simulated times between elongation stages will be short and vice versa.

Once the boundary conditions have been applied, the equilibrium of the system is deduced by means of an algorithm that balances the distribution of ED throughout the domain. The resolution of this balance has been implemented using the iterative algorithm of finite differences presented in Eq. (1). This equation calculates the ED transmitted to each element u_i,j of the matrix in the calculation step (s + 1), based on the values of the ED for each point and their closest neighbours in the previous step, s. If any of the neighbours does not exist, as occurs in the nodes bordering free surfaces or cracks, the amount of equivalent displacement transferred to the considered neighbour is zero. Thus, the equilibrium occurs when the ED gradient is minimal at each point.

$u_{i,j}^{s+1}=u_{i,j}^s+\beta \left(u_{i+1,j}^s+u_{i-1,j}^s+u_{i,j+1}^s+u_{i,j-1}^s-4u_{i,j}^s\right)$  (1)

The β parameter is the ED distribution coefficient and determines the magnitude of the ED transference in each calculation step. A value of β=0.20 has been used. This value is the biggest one that results on a stable behaviour of Eq. (1) and results on the fastest reaching of the equilibrium condition of the whole domain. As a convergence criterion, it is considered that equilibrium is reached when the ED transfer between any of the elements in the matrix is lower than a predetermined value of 10^-6 that has been found to be offer a proper trade-off between calculation time and convergence (Artola, 2018Artola, G. (2018). Susceptibilidad a la fragilización por hidrógeno de aceros de alta resistencia: comportamiento en ambientes marinos y modelización de patrones de agrietamiento. PhD. Thesis, Tecnun-Universidad de Navarra, Donostia-San Sebastián. http://hdl.handle.net/10171/56315.). This type of criterion is required for the method as the equilibrium calculation is not deterministic but asymptotic.

During iterative calculations towards the equilibrium, displacement gradients are generated in the horizontal and vertical directions due to geometrical features in the simulated specimen, whose horizontal and vertical gradient matrices, H_MxN and V_MxN respectively, are calculated according to the Eq. (2) and (3).

$h_{i,j}=\frac{u_{i+1,j}-u_{i-1,j}}{2}$  (2)

$v_{i,j}=\frac{u_{i,j+1}-u_{i,j-1}}{2}$  (3)

From the horizontal and vertical gradient matrices, H and V, operating term by term as indicated in Eq. (4), a scalar that represents the strain state in each point of the domain is calculated.

${def}_{i,j}=\sqrt{v_{i,j}^2+h_{i,j}^2}$  (4)

The DEF_MxN matrix, whose elements are def_i,j , is interpreted as the maximum deformation supported at each point, without considering its direction. From the DEF_MxN, and assuming a linear elastic material from here on, an approximate elastic stress existing at each point of the simulated system can be deduced. In the 2D case studied in this work it was assumed a linear relation between stress and strain. The material Young modulus used is 210 GN·m^-1. As proposed model is 2D, force/length was used instead of force/area that would correspond to the full 3D real system. In this case the strain is applied in the vertical direction. Due to the Poisson’s ratio it is expected to have a horizontal deformation in our virtual sample, but as this Poisson’s related deformation does not produce any stress it was not considered. To know the overall stress applied in a simulated specimen, the stress of all the elements along the top horizontal row of the matrix has been used.

In other words, all the displacement is initially forced on the top row of material and the system is left to find its own equilibrium by distributing (balancing) the displacement among all the nodes in the domain of interest and the gradient of displacement at each point is used to deduce the strain distribution. In the presence of geometrical features that alter the ED distribution, such as in Fig. 2, strain concentrations arise. Assuming a linear relation between strains and stresses it is possible to deduce that the proposed algorithm is capable of reproducing stress concentration values for Fig. 2 configuration, finite-width thin element with a circular hole, that differ 10% with analytically calculated solutions (Peterson, 1974Peterson, R.E. (1974). Stress Concentration Factor. John Wiley and Sons, NY, USA.).

Figure 2.  Example of the displacement model applied to a finite width specimen with a hole in the middle, showing the displacement (left) and approximate strain (right) distributions

A resistance value, r_i,j , is assigned to each element of the domain to determine if a crack grows. This resistance is written in terms of the maximum strain value, def_i,j , each point can stand before breaking. At those points where the resistance value is less than the calculated stress value (always in the condition of a balanced standard and assuming a linear elastic behaviour), the main criterion for crack growth is met. However, the incorporation of some complementary criteria in the algorithm is required, both for the cracking nucleation process and for its growth. These criteria are based on the consideration of a Moore-type environment, i.e, considering the neighbouring elements of face and vertex:

The extra requisites for crack nucleation have been:

• Two contiguous elements (horizontal-vertical or diagonal) must simultaneously overcome the breaking stress.

• Only one of these elements must be in contact with a free surface, that is, the crack must be nucleated in the surface of the material.

• The number of solid neighbors of the one in contact with the surface must be 4 or 5.

• The extra requisites for crack growth after nucleation are:

• Cracks grow cell by cell and strain distribution is reassessed in each growth to determine the length the crack grows due to a given displacement increment during SSRT.

• The number of solid neighbors of the element that has exceeded the critical stress condition must be 7.

• Cracks cannot grow by turning 90 degrees or more.

2.3 Integrated embrittlement model

 ⌅

The introduction of the HE into the model has been carried out by modifying the resistance matrix R_MxN as a function of the concentration of diffusible hydrogen. The interaction effects of the hydrogen with the structural features, such as grain boundaries and dislocations, have been put together in a resistance matrix modification rule. The resistance reduction of each cell is carried out linearly between the following extreme conditions:

• The resistance of the material is considered to fall linearly from 0% to 50% of its original resistance depending on the percentage of the diffusible hydrogen in solution, 0% and saturation condition, respectively. The maximum percent­age of reduction has been taken arbitrarily from the maximum elongation reduction observed in the experimental work used for comparison (Artola et al., 2018Artola, G., Arredondo, A., Fernández-Calvo, A., Aldazabal, J. (2018). Hydrogen embrittlement susceptibility of R4 and R5 high strength mooring steels in cool and warm seawater. Metals 8 (9), 700. https://doi.org/10.3390/met8090700.).

• When all the hydrogen in solid solution is trapped hydrogen, the original strength of the material is considered to remain intact.

The initial resistance has been defined at each point by two components, one fixed and one random. The 2D simulations presented in the Results section of this work have been performed with an average resistance base value of 1 GN·m^-1 (fixed part) and an extra resistance distributed as random uniform component of 10% of the base value. The resulting resistance of cells is in average 1.05 ± 0.05 GN·m^-1. These values mean that base material cells will crack when the calculated stress is around 1 GN·m^-1. As the aim of the model is to give qualitative results for understanding the underlaying physical processes these values simple, but relatively close to real material, were arbitrary chosen.

The correction of the resistance values, due to the diffusible hydrogen, has been performed after each diffusion calculation, so that each simulation step of the cracking patterns for SSRT comprises the following steps:

1. Application of the ED increase in the upper boundary of the domain.

2. Calculation of the ED equilibrium condition and the corresponding approximate strain matrix.

3. Application of the diffusion model on the simulated domain.

4. Correction of the resistance matrix as a function of diffusible hydrogen.

5. Verification, element by element, of the crack growth conditions using the resistance matrix of the previous stage. It is this stage where the cracks actually can grow if all the conditions for it to happen are met.

This work sequence allows defining in each iteration all the property maps that define the model, as shown in Fig. 3.

Figure 3.  Maps calculated by the proposed model: Displacements (U), stresses (DEF), trapped H and diffusible H (C), resistance (R). Domain size 100x400 cells (equivalent to 3.3x13.2 mm).

In order to assess the feasibility of the hydrogen trapping influencing the cracking pattern during SSRT on specimens submitted to HE, the model described above has been used to perform four simulation scenarios. Each of these scenarios corresponded to different capacities of hydrogen trapping: 0%, 25%, 50% and 75% of the total hydrogen solubility. This “hydrogen trapping” concept is an instrumental parameter that stands for the percentage of the total hydrogen concentration that can be trapped in an element. The aim of this parameter is not mimicking the actual trapping mechanism but resembling different trapping trends. As the strength reduction rule due to diffusible hydrogen has been the same for all the four scenarios, any difference in cracking patterns between them is interpreted to be due to apparent diffusion coefficient variations, not to an intrinsic resistance increase of the matrix to diffusible hydrogen. The simulated geometries corresponded to rectangular geometry of 250x800 elements, in which the diffusion of hydrogen from the left and right edges has been allowed. For this case, the cell size was fixed at 33 µm. A fixed hydrogen concentration of 1 has been imposed on both lateral surfaces.

To simulate the increasing stress on the material due to SSRT, a displacement increase has been applied on the top row of the U_MxN matrix at a speed of 4×10^-4 mm per simulation step, which is equivalent to a deformation increase of 5×10^-7 per step. Assuming that the simulated time during loading steps is 25 ms, the simulated strain rate is 2.0×10^-5 s^-1.

Hydrogen diffusion has been simulated using a diffusion coefficient for H in pure recrystallized iron of 1.09×10^-8 m²s^-1 and applying the diffusion algorithm 5 times with a value of α = 0.05 for each step of displacement increase. With all these parameters fixed, it is possible to deduce, using Eq. (6), that the simulated time after applying 5 times de diffusion algorithm is 27.3 ms. To assess the repeatability of the results, 4 different initial random resistance distributions have been simulated for each hydrogen trapping condition. The hydrogen inflow has been set to be continuous across the free surfaces, both the original surface and new crack surfaces created during the simulation. For simulating the hydrogen diffusion through cracks, when a solid cell breaks it is not part of the solid anymore. Its hydrogen concentration is changed to get the same constant value as the rest of surface cells, i.e. to 1. The supplementary material of this paper includes an example of evolution of all maps calculated with the proposed model during the virtual text.

3. RESULTS

Figure 4 presents the results corresponding to three simulated stages of crack growth for specimens with hydrogen trapping capacities of 0%, 25%, 50% and 75% of the total hydrogen solubility. These percentages were chosen as rough representatives of the following conditions and are blind to more precise trapping description approaches such as (McNabb and Foster, 1963McNabb, A., Foster, P.K. (1963). A new analysis of the diffusion of hydrogen in iron and ferritic steels. Institute of Metals Division, AIME 227, 618-627.; Oriani, 1970Oriani, R.A. (1970). The diffusion and trapping of hydrogen in steel. Acta Metall. 18 (1), 147-57. https://doi.org/10.1016/0001-6160(70)90078-7.), meaning they are a parameter aimed to combine, at least, the added effect of the density and energy of each trap type in each material:

Each of these stages matches a stress level beyond the onset of the crack initiation: 100%, 80% and 50% of the maximum stress measured in the top row of the simulated volume during the simulation. It is observed that in cases where low trapping capacity is imposed, 0% and 25%, the nucleation of cracks is higher from the beginning and the hydrogen readily diffuses through the walls of the growing cracks. The other two cases, 50% and 75%, show a lower crack density and as less pronounced diffusion of hydrogen through the fresh fractured surfaces.

According to these results, the simulation supports the hypothesis of a higher trapping capacity being a feasible explanation for the experimental observations of a lower number of cracks in SSRT. This is supported by the fact that the results on Sample B for medium mild cathodic protection (-850 mV) showed a lower HE susceptibility than Sample A specimens, what is known to be associated to a higher hydrogen trapping capacity (Yamasaki and Bhadeshia, 2006Yamasaki, S., Bhadeshia, H.K.D.H. (2006). M₄C₃ precipitation in Fe-C-Mo-V steels and relationship to hydrogen trapping. Proceedings of the Royal Society A 462 (2072), 43-50. https://doi.org/10.1098/rspa.2006.1688.; Fielding et al., 2014Fielding, L., Song, E., Han, D., Bhadeshia, H.K.D.H., Suh, D.-W. (2014). Hydrogen diffusion and the percolation of austenite in nanostructured bainitic steel. Proceedings of the Royal Society A 470 (2168), 20140108-1-17. https://doi.org/10.1098/rspa.2014.0108.).

The repetition of the simulation scenarios employing different random strength distributions allowed to discard these observations to be caused by a probabilistic coincidence. Figure 5 shows the number of cracks that appeared for each trapped hydrogen value, including the dispersion produced by randomized strength. The points on the graph correspond to the means of the 4 simulations performed for each condition and the error bars correspond to the dispersion of these values.

The normalized distance between cracks, that is, the distance between the simulated cracks divided by the width of the simulated geometry was shown in Fig. 6. These numerical results are aligned with the interpretation of the influence of diffusion conditions and hydrogen solubility on the cracking pattern. The patterns with high crack density respond to models with low hydrogen trapping capacity, while the patterns with low crack density correspond to models with very high hydrogen trapping material. It must be noted that, being the instantaneous trapping criteria possibly underestimating the hydrogen diffusion front (Charles et al., 2019Charles, Y., Gaspérini, M., Fagnon, N., Ardon, K., Duhamel, A. (2019). Finite element simulation of hydrogen transport during plastic bulging of iron submitted to gaseous hydrogen pressure. Eng. Fract. Mech. 218, 106580. https://doi.org/10.1016/j.engfracmech.2019.106580.), the trends in the effect of trapping capacity should remain analogous.

In order to finish the discussion on the simulation results, it is very illustrative to build the Stress-Strain curve of each crack growth stage. For this purpose, Fig. 7 has been plotted, showing this stress-strain evolution for the four conditions indicated above. This graph shows that the point corresponding to the maximum strain shifts to the left and the maximum stress falls, as the number of traps is reduced. It is worth underlining that these results are only qualitatively in agreement with the working hypothesis arising from the experimental observations in (Artola et al., 2018Artola, G., Arredondo, A., Fernández-Calvo, A., Aldazabal, J. (2018). Hydrogen embrittlement susceptibility of R4 and R5 high strength mooring steels in cool and warm seawater. Metals 8 (9), 700. https://doi.org/10.3390/met8090700.), as no experimental trapping capacity measurements were available for this work.

4. CONCLUSIONS

Regarding the feasibility of the hydrogen trapping capacity to be a driver for cracking pattern variations between two samples corresponding to the same steel strength grade R4, the simulations agree qualitatively with the experimental observations: