MOLECULAR DYNAMICS SIMULATIONS OF HYDROGEN EMBRITTLEMENT: PRELIMINARY RESULTS. Javier Sanchez1, Pedro de Andres2, Carmen Andrade1, José Fullea1 1 Instituto de Ciencias de la Construcción Eduardo Torroja (CSIC). C/ Serrano Galvache, 4. 28033 Madrid, España. E-mail: javiersm@ietcc.csic.es, andrade@ietcc.csic.es, fullea@ietcc.csic.es 2 Instituto de Ciencias de Materiales de Madrid (CSIC). E-28049 Cantoblanco, Madrid, España. E-mail: pedrodeandres@icmm.csic.es ABSTRACT Hydrogen embrittlement is believed to be one of the main reasons for cracking of structures under stress. High strength steels in these structures often include a ferritic core made of alpha-iron (body centered cubic lattice). Previous work [1] was concerned with the interaction of atomic hydrogen with iron using first principles calculations. We studied the effect of interstitial hydrogen in the iron lattice and the stress induced by the interstitial hydrogen in the host lattice. In this paper we study the dynamical behaviour of hydrogen inside the iron lattice. Using ab-initio Molecular Dynamics we obtain hydrogen diffusion paths and by taking statistical averages we extract diffusion coefficients from Einstein’s equation. Depending on temperature, the diffusion path involve going through tetrahedral or octahedral sites. Simulations where a number of hydrogens occasionally coincide in one unit cell have been performed to elucidate the effect of interactions between hydrogens. . KEY WORDS: Hydrogen Embrittlement, Molecular Dynamics, alpha-iron. 1. INTRODUCTION Los aceros de alta resistencia empleados en estructuras civiles o estructuras singulares se componen de una matriz ferrítica, o lo que es lo mismo de una estructura de hierro cúbica centrada en el cuerpo (BCC). Uno de los más empleados a nivel mundial es el acero trefilado de composición eutectoide. Diversos estudios muestran el crecimiento de fisuras para valores de tensión inferiores a su correspondiente tenacidad de fractura bajo la acción de un medio agresivo. Este fenómeno, conocido como Corrosión Bajo Tensión (CBT), es uno de los problemas que afectan a estos aceros. Existe un convencimiento general de que el hidrógeno juega un papel importante en este proceso, de hecho se supone que la fragilización por hidrógeno es una de las causas más frecuentes de fallo en estructuras sometidas a esfuerzos mecánicos [1]. Teniendo en cuenta esta hipótesis se desarrolló un ensayo basado en una disolución de tiocianato amónico para determinar la susceptibilidad de los aceros de alta resistencia a la Fragilización por Hidrógeno. Este ensayo se denominó FIP-78, y fue propuesto por la Federación Internacional de Pretensado. Existe evidencia de la variación de parámetros mecánicos [2], pero la fragilización por hidrógeno no está completamente explicada desde un punto de vista teórico. En este sentido, varios modelos tratan de explicar la propagación de la fisura por la presencia en el metal de átomos de hidrógeno. Generalmente se asume que el hidrógeno se genera electroquímicamente en la superficie del material y difunde hasta la zona en proceso de fractura. Para explicar el proceso por el cual el hidrógeno fragiliza el material existen varias teorías: • Cambio estructural o de fase producido por el hidrógeno [3, 4, 5]. • Plastificación producida por el Hidrógeno o hydrogenenhanced localizad plasticity (HELP) [6, 7]. • Reducción de la energía cohesiva por el efecto del hidrógeno [8, 9, 10]. Para entender mejor estas teorías, hemos realizado cálculos atomísticos de primeros principios [11] con el objetivo de determinar la posición preferente del hidrógeno intersticial y las tensiones que genera en la red bcc del hierro en función de la densidad de hidrógeno en la red de hierro. Dado que estos cálculos son computacionalmente muy costosos cuando incluyen muchos átomos en la celda unidad, hemos extendido estos resultados por medio de cálculos de Elementos Finitos. Boltzmann y T es la temperatura absoluta. En este trabajo se presentan los resultados preliminares obtenidos a través de cálculos de Dinámica Molecular Ab-Initio. Esta técnica es complementaria a los resultados obtenidos anteriormente y permite observar los procesos dinámicos que tienen lugar dentro de la red de Fe con diferentes concentraciones de H. 3. RESULTADOS Los cálculos de Dinámica Molecular Ab-Initio se han llevado a cabo con el programa CASTEP [12,13]. Se han realizado en una supercelda 2x2x2 incluyendo 16 átomos de Fe y en la cual se han introducido diferente numero de átomos de H según el caso. De esta forma es posible estudiar el efecto de la concentración de H. Se ha considerado la aproximación de Born-Oppenheimer que aplica los principios de la mecánica clásica a los iones, los cuales son objetos que se mueven en potencial creado por los electrones, que se tratan como objetos quánticos y obedecen a la ecuación de Schroedinger. Las funciones de onda de Kohn-Sham se han desarrollado en una base de ondas planas. La precisión de los cálculos viene determinada básicamente por dos parámetros: (i) la máxima energía de corte (“cut-off”) que en este caso es de 375eV y (ii) el número de puntos usados en la zona de Brillouin para muestrear las funciones de onda en espacio reciproco (“puntos k”), que en este caso se ha considerando una malla 4x4x4 tipo Monkhorst-Pack [MP]. Se han utilizado pseudopotenciales ultra-suaves [14] y la aproximación de gradientes generalizados de Perdew, Burkeand Ernzerhof [15]. En el estudio de la difusión de H dentro de la red de Fe hemos realizado simulaciones previas en el conjunto micro-canónico (la energía total es una constante del movimiento y el volumen de la celda es constante). Se han realizado simulaciones de 1-2ps con pasos de tiempo de 0.5-1.0fs, en los cuales la energía total se conserva con un error de 0.01%. De los cálculos se han excluido los resultados de los primeros 100-200fs en los cuales el sistema se está equilibrando. Los coeficientes de difusión se han calculado asumiendo que los átomos intersticiales de H se mueven de forma aleatoria (random-walk). Se han trazado las trayectorias de los átomos de H durante el tiempo de la simulación y se han calculado los desplazamientos para obtener el coeficiente de difusión (D) de acuerdo a la ecuación de Einstein en tres dimensiones: <|r(t) - r(0)|2> = 6Dt. donde: r es la posición del átomo y t es el tiempo. (1) Las barreras de difusión se han obtenido a través de la ecuación de Arrhenius: (2) D = D0*exp[-Ea/kB*T] donde D0 es el factor pre-exponencial, Ea es la energía de activación o barrera de difusión, kB es la constante de Fe H Figura 1. Trayectoria del H (puntos negros) dentro de la red de Fe (circulo rojo). T=300K. En la Figura 2 se muestra el camino recorrido por el H a lo largo del tiempo. A través de la ecuación de Einstein (ecuación 2) es posible estimar el coeficiente de difusión para el H. 1 <|r(t) - r(0)|2> 2. METODOLOGÍA A continuación se presentan los resultados obtenidos de los cálculos de Dinámica Molecular para diferentes concentraciones de H. En primer lugar se muestran los resultados obtenidos para el caso de 1 átomo de H en la supercelda 2x2x2 de Fe, es decir, una relación H/Fe=1/16. En la Figura 1 se muestra la trayectoria del H dentro de la red de Fe para una temperatura de 300K. Para esta concentración y esta temperatura el camino de difusión del H es a través de los huecos tetraédricos, evitando los huecos octaédricos. difs( x) 0.5 FF ( x) 0 0.5 0 500 1000 1500 Timex (fs) Figure 2. Random-walk para un H a 300K. En la Figura 3 se muestra la trayectoria del H dentro de la red de Fe para una temperatura de 700K. Al aumentar la temperatura aumenta la movilidad del átomo de H de acuerdo con la ley de Boltzmann (Eq. 2). Al mismo tiempo, cambia el camino de difusión: a 300 K el átomo de H reside aproximadamente un 0.65% del tiempo cerca del hueco octaédrico, mientras que a 700 K permanece algo mas del 6%. Este efecto se puede entender fácilmente calculando la función de partición de un sistema con dos niveles y al igual que el aumento de movilidad solo depende de la relación entre la energía cinética, kT, y la diferencia de energía entre los dos sitios (o la barrera en el caso de la difusión). 1.E-06 D (m2/s) 1.E-07 1.E-08 1.E-09 0.001 0.0015 0.002 0.0025 0.003 0.0035 (1/T) Figure 5. Arrhenius plot. Ea= 145meV Figure 3. Trayectoria del H (puntos negros) dentro de la red de Fe (circulo rojo). T=700K. En la Figura 4 se muestra el camino recorrido por el H a lo largo del tiempo a 700 K. Los cálculos indican un incremento del coeficiente de difusión desde 5.6*109 2 m /s a 300K hasta 2.3*10-7m2/s a 700K. En la Figura 6 se muestran los coeficientes de difusión para el caso de una concentración de H/Fe = 8/16. En este caso, para cada cálculo de Dinámica Molecular se obtienen 8 valores del coeficiente de difusión, uno por cada átomo de H: los resultados presentados corresponden al promedio de todos ellos. De acuerdo con la ecuación de Arrhenius (ecuación 2) el valor de la Ea obtenido es de 47 meV, concluyendo que la barrera de difusión es afectada por el número de intersticiales difundiendo simultáneamente. 2 1.5 difs( x) FF ( x) 1 0.5 0 0.5 Hasta ahora se ha estudiado el efecto de la temperatura sobre la difusión del H en la red de Fe. La segunda variable de interés es la concentración de H dentro de la red de Fe. Se han repetido los cálculos para diferentes concentraciones: H/Fe = 2/16, 4/16 y 8/16, y para diferentes temperaturas: T = 300, 400, 500, 600 y 700K 0 500 1000 1500 x Figure 4. Random-walk para un H a 700K. En la Figura 5 se han representado los coeficientes de difusión obtenidos de los cálculos de Dinámica Molecular para diferentes temperaturas y una concentración H/Fe=1/16. Teniendo en cuenta la ecuación de Arrhenius (ecuación 2) es posible obtener la Energía de Activación (Ea) correspondiente a la difusión del H dentro de la red de Fe. En este caso (baja concentración de H) se obtiene una Ea= 145 meV. Esta energía de activación está próxima a la calculada por Primeros Principios en trabajos anteriores [11,16]. En nuestro trabajo previo hemos concluido que la estabilidad de los sitios de alta simetría puede depender de la densidad de H y las restricciones de simetría que se utilicen como condiciones de contorno [11,16]. Para baja concentración de H el sitio mas estable es el tetraédrico mientras que para una concentración de ½ el hueco octaédrico se convierte en mas estable y se favorece una distorsión tetragonal de la simetría cúbica. Este aspecto se ha confirmado con los cálculos de Dinámica Molecular, donde se observan dos hechos relevantes: (i) se produce un cambio en el camino de difusión para baja concentración de H (de tetraédrico a tetraédrico), a un camino de difusión para alta concentración de H que va de octaédrico a octaédrico pasando por el tetraédrico. En segundo lugar, (ii) se produce una disminución de la barrera de difusión (Ea) o, lo que es lo mismo, se produce un aumento del coeficiente de difusión. 1.E-06 D (m2/s) 1.E-07 1.E-08 1.E-09 0.001 0.0015 0.002 0.0025 0.003 0.0035 (1/T) Figure 8. Vibración del átomo de Fe alrededor de su posición de equilibrio. Figure 6. Arrhenius plot. E=47meV Ocupación hueco octaédrico (%) En la Figura 7 se muestra cómo aumenta la ocupación de los huecos octaédricos al aumentar la concentración de H. En este caso se muestran los resultados correspondientes a una temperatura de 500K. El efecto de la concentración de H sobre la ocupación de los huecos octaédricos es superior al efecto de la temperatura porque la variación en la diferencia de energía de los dos sitios de simetría es más importante que la variación de energía térmica al pasar de 300 a 700 K. 60 50 40 30 20 10 0 0 2 4 6 8 10 Numero de H en la supercelda 2x2x2 Figure 7. Ocupación de los huecos octaédricos frente a la concentración de H para T=500K. 3.1. Comportamiento de los átomos de Fe. Otro aspecto a analizar es el comportamiento de los átomos de Fe en presencia de H. Los átomos de Fe se mueven alrededor de su posición de equilibrio. El desplazamiento cuadrático medio <u2> de los átomos de Fe se puede calcular en un modelo isótropo de Debye en función de un único parámetro ajustable, la temperatura de Debye (Θ), En la Figura 9 se muestra la distribución de desplazamientos ajustada a una distribución normal centrada en la posición de equilibrio y cuya varianza coincide con el desplazamiento cuadrático medio (σ2 ~ <u2>). Como para el caso del coeficiente de difusión, se ha estudiado la influencia de la concentración de H sobre el comportamiento del Fe, ya que el comportamiento frente a la temperatura viene definido por el modelo de Debye-Waller (ecuación 3). Como se puede observar en la Figura 10, el desplazamiento cuadrático medio de los átomos de Fe respecto a su posición de equilibrio aumenta de forma considerable al aumentar la concentración de H intersticial. En este caso se muestran los resultados obtenidos a 700K, pero este comportamiento se repite en todo el rango de temperaturas estudiado. El aumento del desplazamiento <u2> de los átomos de Fe indica el debilitamiento del enlace de los mismos por la presencia del H. Esto está de acuerdo con los diversos resultados experimentales existentes en la literatura sobre la Fagilización por Hidrógeno y ensayos de Corrosión Bajo Tensión [1728], donde el efecto del H se manifiesta con la reducción del tiempo de rotura y la reducción de propiedades mecánicas, como por ejemplo la tenacidad de fractura. 300 10 8 200 6 eventos eventos_sim 4 100 2 (3) En la Figura 8 se muestra un ejemplo del desplazamiento de un átomo de Fe alrededor de su posición de equilibrio durante un cálculo de Dinámica Molecular. 0 0 0.1 0.2 0 0.3 desplazamiento Figure 9. Desplazamiento de los átomos de Fe respecto su posición de equilibrio. 0.08 Materialia 51 (2003), pp. 2717–2730. y = 0.00367x + 0.04851 R2 = 0.98034 0.075 [3] DG Westlake - Trans. Metall. Soc. AIME, 1969. <u> prom Fe (A) 0.07 [4] Nelson H.G. “Film-rupture model of hydrogeninduced, slow crack growth in acicular alpha-beta titanium” Metall. Trans. A - Phys. Metall. Mater. Sci. 7A (5), pp. 621–627, 1976. 0.065 0.06 0.055 0.05 0.045 0.04 0 1 2 3 4 5 6 7 8 9 Num H en la celda 2x2x2 Figure 10. Desplazamiento promedio de los átomos de Fe en función de la concentración de H para T=700K. 4. CONCLUSIONES i. El proceso de difusión del H depende de la temperatura de una forma estándar a través de un factor de Boltzmann, pero también de la concentración de H intersticial. ii. El aumento de la concentración de H produce un cambio en el camino de difusión: para bajas concentraciones el H difunde a través de huecos tetraédricos contiguos, mientras que para alta concentración de H la difusión se produce a través de huecos octaédricos pasando por el tetraédrico que queda entre ellos. iii. El aumento de la concentración de H produce la disminución de la barrera de difusión, o lo que es lo mismo, produce el aumento del coeficiente de difusión. iv. La presencia de los átomos de H dentro de la red de Fe producen un debilitamiento del enlace, lo cual puede llegar a modificar la resistencia mecánica de los aceros de alta resistencia. Estos resultados obtenidos con los cálculos de Dinámica Molecular permiten entender mejor los resultados experimentales disponibles. [5] Oriani R.A. “Hydrogen Effects in High-strength Steels” In: Gangloff, R.P., Ives, M.B. (Eds.), First International Conference on Environment-induced Cracking of Metals, NACE-10. NACE, Houston, TX, pp. 439–447, 1990. [6] Beachem, C.D., 1972. “A new model for hydrogenassisted cracking (hydrogen embrittlement)”. Met. Trans. 3, 437–451. [7] Birnbaum, H.K., Sofronis, P. “Hydrogen-enhanced localized plasticity – a mechanism for hydrogen related fracture”. Mater. Sci. & Eng. A 176, 191– 202., 1994. [8] Troiano A R 1960 Trans. ASM 52 54. [9] Oriani R.A. “Mechanistic theory of hydrogen embrittlement of steels. Ber. Bunsenges” Phys. Chem. 76 (8), pp. 848–857, 1972. [10] Oriani, RA; Josephic, PH, "Equilibrium and kinetic studies of the hydrogen-assisted cracking of steel" Acta Metall. 1977. pp. 979-988. [11] J. Sanchez, J. Fullea, C. Andrade, and P. de Andres, Phys. Rev. B 78, 014113 (2008). [12] S. Clark, M. D. Segall, C. Pickard, P. Hasnip, M. J. Probert, K. Refson, and M. C. Payne, Z. fuer Kristal-lographie 220, 567 (2005). [13] (CASTEP 4.4); http://www.accelrys.com. [14] D. Vanderbilt, Phys. Rev. B 41, 7892, 1990. ACKNOWLEDGEMENTS This work has been financed by the Spanish CICYT (MAT2008-1497), and MEC (CONSOLIDERS CSD2007-41 ”NANOSELECT” and ”SEDUREC”). REFERENCES [1] M. Elices “Influence of residual stresses in the performance of cold-drawn pearlitic wires” Journal of Materials Science, Volume 39, Number 12, 2004, pp. 3889 – 3899. [2] Y. Liang, P. Sofronis and N. Aravas “On the effect of hydrogen on plastic instabilities in metals” Acta [15] J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996). [16] J. Sánchez, P. de Andrés, J. Fullea, C .Andrade. ”Aproximación por simulación ab-initio a la fragilización por hidrógeno en una red de hierro bbc”. Anales de Mecánica de la Fractura, Vol. 24, 387-392, 2007. [17] Toribio J., Elices M. “Nuevas Aportaciones al Ensayo FIP de Fragilización por Hidrógeno en Tiocianato Amónico” Hormigón y acero, Vol. 27 (168), pp. 121-130, 1988. [18] Pakins R. N., Zhou S. “The Stress Corrosion Cracking of C-Mn Steel in CO2-HCO3– - CO32Solutions I: Stress Corrosion Data” Corr. Sci., Vol. 39, Nº 1, pp. 159-173, 1997. [19] Caballero L., Elices M. “Influencia de la velocidad de deformación en la propagación de fisuras por corrosión bajo tensión.” Revista Iberoamericana de Corrosión y Protección, 17 (1), pp. 15-22, 1986. [20] Caballero L., Elices M. “Un método para la medida de la cinética de las fisuras de corrosión bajo tensión en ensayos a velocidad de deformación constante.” Revista Iberoamericana de Corrosión y Protección, 17 (1), pp. 43-48, 1986. [21] Acha-Hurtado M. “Corrosión Bajo Tensión de Alambres de Acero Pretensado en Medios Neutros con HCO3– y Alcalinos con SO4=.” Instituto de Ciencias de la Construcción “Eduardo Torroja” C.S.I.C. Madrid, 1993. PhD Thesis. [22] Alonso M. C., Andrade C., Procter R. P. M., Saenz de Santa María M. “Susceptibilidad a la Corrosión Bajo Tensión del Acero Pretensado en Disoluciones de NaHCO3” Hormigón y Acero, Nº 166, pp. 121126, 1988. [23] Lancha A. M. “Influencia del Trefilado en la Corrosión Bajo Tensión de Aceros Eutectoides” Universidad Complutense de Madrid. Facultad de Ciencias Químicas. 1987. [24] Valiente A., Elices M. “Premature Failure of Prestressed Steel Bars” Engineering Failure Analysis, Vol. 5, nº 3, pp. 219-227, 1998. [25] Sanchez J., Fullea J., Andrade C., Alonso C. “Stress corrosion cracking mechanism of prestressing steels in bicarbonate solutions” Corr. Sci. 49(11), 4069-4080, 2007. [26] Sanchez J., Fullea J., Andrade C. “Fracture toughness variation induced by stress corrosion cracking of prestressing steels” Materials and Corrosion, 59(2), 2008. [27] Sanchez J., Fullea J., Andrade C. “Reasons for Crack Arrest in Stress Corrosion Cracking Tests. Crack Propagation Rate in High Strength Steels” Corrosion, 65(6), 2009. [28] Sanchez J., PhD, Thesis. Eduardo Torroja Institute, IETcc-CSIC, Spain, 2007.