Anales AFA Vol. 32 Nro. 4 (Enero 2022 – Abril 2022) 112-119
https://doi.org/10.31527/analesafa.2021.32.4.112
Materia Condensada
ESTUDIO TEÓRICO Y NUMÉRICO DE LA DIFUSIÓN DE HIDRÓGENO EN
ALEACIONES BASE Zr-Nb
THEORETICAL AND NUMERICAL STUDY OF HYDROGEN DIFFUSION IN Zr-Nb
ALLOYS
G. Beltramo1, V. P. Ramunni2,3 y R. C. Pasianot*1,2,3
1 Instituto Sabato, UNSAM/CNEA Av. Gral. Paz 1499, (1650) San Martín, Argentina.
2 Consejo Nacional de Investigaciones Científicas y Técnicas – CONICET, Godoy Cruz 2390 (C1425FQD),
Argentina.
3 Gcia. Materiales, Comisión Nacional de Energía Atómica, Av. Gral. Paz 1499, (1650) San Martín, Argentina.
Autor para correspondencia: email: pasianot@cnea.gov.ar
ISSN 1850-1168 (online)
Recibido: 02/06/2021 Aceptado: 15/08/2021
Resumen:
En el marco de los problemas de fisuración asistida por hidruros en los tubos de presión de los reactores nucleares
tipo CANDU, estudiamos la influencia de la fase continua Zr-β (∼ 20% Nb - ∼ 80% Zr) sobre el coeficiente de
difusión del hidrógeno, DH, en las aleaciones bifásicas Zr-α / Zr-β. Proponemos un modelo fenomenológico de la
difusión en tubos que constituye una mejora respecto de otros disponibles en literatura. También estudiamos la
influencia del contenido de Nb sobre el DH en la fase Zr-β (cúbica) empleando la teoría del estado de transición,
provista de parámetros ab-initio calculados con el código SIESTA. En particular, consideramos 9 aleaciones
ordenadas con distintas concentraciones de Nb, y determinamos energías de activación efectivas mediante ajustes
de Arrhenius para cada aleación. Encontramos que la concentració de Nb varía la energía de activación de forma
no monótona, pasando por un máximo para Zr-50%Nb. Finalmente, observamos que la mayor difusividad en la
dirección longitudinal (o axial) del tubo vs. la radial, predicha y medida, es consistente con la microestructura
texturada del material. Concluimos que la pérdida de continuidad de las láminas de Zr-β presentes en la
microestructura de los tubos de presión, es consistente con la disminución de DH a medida que trascurre el tiempo a
una dada temperatura.
Palabras clave: Zr-Nb, Difusión de H, Ab-Initio, fisuración asistida por hidruros.
Abstract:
Within the framework of the delayed hydride cracking phenomena reported for the pressure tubes of CANDU-type
nuclear reactors, we study the influence of the Zr-β (∼20%Nb - ∼80%Zr) continuous phase on the hydrogen
diffusion coefficient, DH, through the Zr-α / Zr-β biphasic alloy.We propose an improved phenomenological model
for DH with respect to those found in the literature. Furthermore, we study the influence of the Nb content on DH
for the cubic phase Zr-β, employing the transition state theory furnished with ab-initio parameters provided by the
SIESTA code. In particular, 9 ordered alloys are considered with different Nb content and effective activation
energies are computed by Arrhenius fits for each alloy. We find that activation energies vary in a non-monotonic
way as Nb content increases, reaching a maximum value at about Zr-50%Nb. Finally, we observe that the predicted
and measured larger diffusivity along the tube axis vs. the radial direction, is consistent with the material texture.
Moreover, we conclude that the loss f continuity of the Zr-β sheets present in the tube microstructure, is consistent
with the decrease of DH in time at a given temperature.
Keywords: Zr-Nb, H Diffusion, Ab-Initio, delayed hydride cracking.
I. INTRODUCCIÓN
Las aleaciones base zirconio son utilizadas en los reactores nucleares debido a su baja absorción de neutrones y
buena respuesta a la corrosión a alta temperatura,1 en particular, la aleación Zr-2,5%Nb es utilizada en la
fabricación de los tubos de presión de las centrales nucleares tipo CANDU (CANadá Deuterio Uranio).2 En éstas,
el hidrógeno altera la mayorıa de las propiedades mecánicas del material y puede llevar a la rotura catastrófica del
componente a través del proceso conocido como fisuración asistida por hidruros (DHC). Teniendo en cuenta que el
hidrógeno puede estar presente en las aleaciones como una impureza remanente del proceso de fabricación y
además puede ingresar durante la vida en servicio del material, es de suma importancia evaluar la respuesta de las
aleaciones frente al hidrógeno.3 Así, a fin de evaluar el daño producido al componente y pronosticar su vida útil en
servicio, resulta fundamental poder caracterizar la difusión del H en el material. En estas condiciones, los tubos son
sometidos a temperatura e irradiación, que modifican sus propiedades mecánicas, microestructura, solubilidad y
coeficiente de difusión del H, lo cual afecta directamente la DHC.
La aleación Zr-2.5Nb es bifásica con granos de fase Zr-
rodeados de láminas de Zr-
.4 La Fig. 1 muestra
micrografías del material de los tubos de presión, obtenidas por microscopía electrónica de transmisión (TEM). Las
láminas originales de Zr-
son metaestables en el rango de temperaturas de funcionamiento del reactor.
FIG. 1: Micrografías TEM del material de tubos de presión tomadas de (a) ref. 4 y (b) ref. 1. La flecha señala
pequeños precipitados de Zr-ω en una matriz de Zr-β enriquecida.
Por esta razón se degradan con el tiempo perdiendo su continuidad, e interrumpiendo ası los caminos
(supuestamente) rápidos para la difusión del hidrógeno en la dirección longitudinal del tubo, al mismo tiempo que
se modifica su contenido de Nb. Los experimentos realizados en nuestro laboratorio5 plantean dos importantes
interrogantes que pretendemos responder en este trabajo, ellos son:
• La pérdida de continuidad de las láminas originales ¿es efectivamente la responsable de la disminución
observada en el coeficiente de difusión del H a medida que trascurre el tiempo a una dada temperatura?
• La difusión y/o solubilidad del hidrógeno en la aleación ¿se ven afectadas por el contenido de Nb de la fase
Zr-
?
En las secciones que siguen presentamos las metodologías que utilizamos para abordar estas cuestiones y los
resultados que de ellas se derivan. Brevemente, i) desarrollamos un modelo fenomenológico de la difusión en los
tubos de presión, con el fin de evaluar la influencia de la fase continua Zr-
sobre el coeficiente de difusión del
hidrógeno, tanto axial como radial; y ii) empleamos la teoría del estado de transición6 (TST), provista de
parámetros calculados ab initio mediante el código SIESTA,7 para estudiar el coeficiente de difusión en el Zr-
conforme varía la concentración de Nb.
II. METODOLOGÍA
Modelo de la difusión en tubos
Proponemos un modelo de difusión unidimensional, Fig. 2, basado en el de Skinner y Dutton (S&D),1 que consiste
de un arreglo tridimensional formado por un paralelepípedo de fase Zr-
, rodeado en tres de sus caras por finas
láminas,
,LR
y
T
, de una fase bifásica, que llamaremos
a secas, compuesta por 31% de Zr-
enriquecida y
69% de Zr-
según los resultados de S&D.1 Esto último correspondería a la condición del tubo luego del
tratamiento ordinario en autoclave a 400 ºC por 24 h, previo al ingreso en reactor (c.f. Fig. 1). Dichas láminas
cumplen la función de caminos rápidos para la difusión del H, en consonancia con el modelo original de Sawatsky
et al.4
Observamos que en cualquier dirección del tubo, el camino efectivo de difusión es una sucesión de granos
seguidos de finas láminas
, en paralelo con un camino continuo compuesto por láminas
únicamente. Por
ejemplo, en la dirección longitudinal el grano
está en serie con la lámina
L
de espesor
L
, y a su vez este
conjunto está en paralelo con las láminas
R
y
T
de espesores
R
y
T
respectivamente. La Fig. 3 muestra el
análogo eléctrico del caso.
FIG. 2: Bloque de difusión modelo, esquemático.
FIG. 3: Circuito equivalente para la dirección
L
.
Considerando que el coeficiente de difusión desempeña el papel de una conductividad, la resolución del circuito
conduce a la siguiente expresión del coeficiente de difusión efectivo en la dirección longitudinal del tubo,
L
D
,
*
(1 ) ,
1 ( 1)
RT L RT
L
RT
D KD
DK
(1)
donde
()
RT R T
y
*
1
1{1 ( 1)}.
LL L
L
K
D D KD
(2)
En las Ecs. (1) y (2),
,
D
son coeficientes de difusión,
/
X X X
d
es la fracción volumétrica de lámina
X
,
K
es el coeficiente de partición de solubilidades del H en fase
respecto de
. Una expresión análoga se obtiene
para
R
D
sin más que intercambiar
R
con
L
.
Resumimos a continuación los datos del modelo. Para las difusividades ,
0( / )D D exp Q RT
, (
[]D
m
2
/s y
[]Q
KJ/mol):
07
7.73 10D
,
45.348Q
de ref.8 y
07
1.1 10D
,
27.139Q
de ref..1 Para el coeficiente
de partición usamos
2K
.4 Los datos geométricos se consignan en la Tabla 1 (c.f. Fig. 2).
TABLA 1: Parámetros geométricos.
Las dimensiones del grano
respetan proporciones típicas para los tubos.9 Las variables 1 y 2 se fijan buscando el
mejor ajuste con las mediciones de cociente de velocidades de propagación de fisuras (
P
V
) en sentido longitudinal
vs. radial; según Dutton et al.10
P
V
es proporcional a las difusividades. Además imponemos la condición que el
volumen de las láminas rodeando al grano
no debe exceder el 10% del total.11
Coeficientes de difusión ab initio
Para el caso escalar, adecuado a la presente situación de redes cúbicas, el coeficiente de difusión
D
adopta la
expresión
2,
6t
t
R
D lim t
(3)
donde
t
R
es la posición del trazador a tiempo
t
y
indica promedio sobre el ensamble de trayectorias.
Considerando al H como un caminante aleatorio por la red BCC,
R
se compone de saltos elementales entre sitios
tetraédricos (T) primeros vecinos, conectados por barreras de energía. En el marco de la TST6 la frecuencia
de
estos saltos se expresa según,
3
/
1
2
1
( / 2 )
2,
( '/ 2 )
iE kT
i
sinh h kT
kT e
hsinh h kT
(4)
donde
E
es la barrera de energía determinada por la configuración de ensilladura entre sitios T y
i
y
'
i
son las
frecuencias de vibración del H en T (equilibrio) y en el punto silla, respectivamente. Obsérvese que para la silla se
suprime la frecuencia imaginaria asociada y además hemos supuesto, dado la diferencia de masas, que el H vibra
desacoplado de la red. A altas temperaturas el factor pre-exponencial se reduce al clásico cociente de frecuencias.
En este contexto, para el caso de una matriz BCC pura (e.g. H en Nb) la Ec. (3) toma la forma explícita,
2
2,
3
H
Ds
(5)
siendo
s
la longitud del salto. Lamentablemente las expresiones analíticas resultan imposibles ni bien la simetría
se degrada, máxime en aleaciones desordenadas como el Zr-
. En este trabajo enfocamos el estudio sobre nueve
aleaciones con distintos contenidos de Nb. Cada una de estas, si bien ordenada, posee un conjunto de saltos no
equivalentes, razón por la cual atacamos el cálculo de
H
D
mediante simulación numérica de la Ec. (3), empleando
la técnica de Monte Carlo cinético12 (KMC), alimentada con frecuencias del tipo Ec. (4) en el límite clásico
(suficiente para nuestros propósitos). Esto produce gráficos de Arrhenius,
H
logD
vs.
1/T
, ligeramente curvos, que
ajustamos con una sola energía de activación
Q
(tomada como representativa) en un intervalo de temperaturas que
comprende las de operación del reactor.
Calculamos todas las frecuencias y energías necesarias mediante el código SIESTA,7 basado en la teoría del
funcional densidad (DFT), pseudopotenciales y bases numéricas localizadas. Los dos últimos aspectos fueron
desarrollados y optimizados en este trabajo, a fin de lograr una pareja de especies Zr-Nb lo más consistente posible.
La mayoría de los cálculos se realizaron sobre (super) celdas cúbicas BCC de 16 nodos (
222
celdillas
unitarias), con distinto contenido de Nb: 0, 1, 2, 4, 8, 12, 14, 15 y 16 átomos, nueve en total, que en lo sucesivo
notaremos como
( /16)n
. Estas celdas fueron generadas con la máxima simetría posible a fin de minimizar el
número de saltos a considerar. Para c/u de ellas encontramos el parámetro de red óptimo, luego lo fijamos e
introdujimos 1 átomo de H en distintas posiciones, que a su vez optimizamos mediante relajación de coordenadas.
Así determinamos que los sitios tetraédricos son mínimos de la energía, en tanto que los octaédricos (de mayor
energía) resultan sillas de rango 2. En cuanto a las sillas entre dos sitios T, la mayoría puede evaluarse ubicando al
H a mitad de camino sobre el segmento entre ambos y minimizando la energía con la restricción de mantener al H
sobre el plano perpendicular al segmento. Para los pocos casos no simétricos recurrimos al método de “drag”,13 que
consiste en levantar el perfil de energía para posiciones iniciales del H interpoladas entre dos configuraciones T
relajadas, y minimizar la energía anulando la componente de la fuerza sobre el (hiper)segmento de interpolación.
En todas estas simulaciones con SIESTA empleamos los siguientes parámetros principales: i) paso de malla
espacial equivalente a 450 Ry (MeshCutoff), ii) temperatura electrónica de 0.15 eV (distribución de Fermi-Dirac),
iii) malla en espacio recíproco de
555
puntos, iv) fuerzas residuales por debajo de 0.01 eV/Å para relajaciones
convergidas.
III. RESULTADOS Y DISCUSIÓN
Modelo fenomenológico
La Fig. 4 muestras difusividades axiales del hidrógeno para tubos bajo distintas condiciones, ya sea medidas, en
líneas llenas, o predichas por modelos, en líneas de trazos. La curva (a) corresponde a Zr-
,8 (b) y (c)
corresponden a Zr-2.5wt%Nb con tratamientos térmicos de 565 ºC y 400 ºC respectivamente,1 (d) es la medición
de 4 también para Zr-2.5wt%Nb, (e) es la predicción de,4 (f) la de1 y (g) la del presente trabajo. Para esta última, las
variables 1 y 2 de la Tabla 1 se fijaron en 0.08, lo que resulta en una fracción de volumen de 9% para la fase
,
por debajo del máximo permitido de 10%11. El rango de temperaturas de 453 a 543 K es el que en general se utiliza
para los ensayos de DHC.11
FIG. 4:
H
D
axial medido para tubos en distintas condiciones, líneas llenas, y predicciones de modelos, líneas de
trazos. Detalles en texto principal.
Observamos que nuestro modelo concuerda muy bien con los resultados experimentales de Sawatsky et al.4
Por otra parte, en la Fig. 5 mostramos el cociente de
P
V
axial/radial medido por Jovanović et al.,14 curvas (a), (b) y
(c), en material Zr-2.5wt%Nb de tubos con tratamientos térmicos a 400 ºC por 1 h, 24 h y 1000 h, respectivamente.
La curva (d) es lo mismo pero según la norma canadiense CSA 285.15;15 finalmente (e) es nuestra predicción para
el cociente de difusividades axial/radial. Esta última tiene un comportamiento parecido al de la norma canadiense
para el cociente de
P
V
y está dentro de los valores consignados en literatura para esta magnitud.1,8 Ambas difieren
sin embargo de los resultados de Jovanović et al., curvas (a) y (b), aunque también dentro del rango de valores. Al
presente, las razones de esta discrepancia son de carácter sumamente especulativo.
FIG. 5: Cocientes de
P
V
axial/radial (a-d) y nuestra predicción para el cociente de
H
D
axial/radial (e).
Respetando la limitación mencionada respecto de la fracción volumétrica de las fases presentes, es posible ajustar
los espesores
L
,
R
, y
T
de nuestro modelo con el fin de alcanzar otros valores máximos y mínimos de cociente
axial/radial para
H
D
. El rango de valores se muestra de forma esquemática en la Fig. 6. Los valores máximos
fueron alcanzados con
0
L
,
0.08
R
, y
0
T
, obteniendo un volumen para las láminas que rodean al grano
de Zr-
de 7.4%; los valores mínimos fueron logrados para
0.3
L
,
0.08
R
y
0.1
T
, obteniendo una
fracción volumétrica de aproximadamente 10%.
FIG. 6: Espacio de valores accesible para el cociente de
H
D
axial/radial.
Cálculos ab initio preliminares
A fin de validar nuestras predicciones, en primer término realizamos cálculos de parámetros de red y energía de
activación para
H
D
en Nb puro, dado la existencia de resultados experimentales confiables así como de cálculos de
literatura. En la Fig. 7 mostramos la variación del parámetro de red en función del contenido de Nb, predicha a
partir de las nueve aleaciones mencionadas, comparada con otras estimaciones ab initio16 mediante el código
Wien2k,17 y una interpolación de resultados experimentales.18,19 Nuestros cálculos emplearon celdas cúbicas de 16
nodos, los de ref.16 celdas no cúbicas de 32 nodos, y los experimentos corresponden a aleaciones desordenadas
retenidas desde temperaturas y concentraciones de Nb que estabilizan la fase BCC (
700
K y
10%‰
resp.).
Ambos cálculos muestran un razonable acuerdo con los experimentos, reproduciendo en particular el
comportamiento lineal tipo ley de Vegard, si bien el nuestro con una subestimación un poco mayor.
FIG. 7: Parámetro de red vs. contenido de Nb; (a) ajuste a datos experimentales,18,19 (b) cálculo de ref.,16 (c)
cálculo presente.
En la Tabla 2 presentamos energías de activación
A
Q
para las especies puras, calculadas y medida en Nb. Nuestros
valores marcados con (*) representan una extrapolación a celda de tamaño infinito, basada en teoría elástica y
siguiendo metodologías de literatura.22,23 Su coincidencia para los dos tamaños que probamos, 16 y 54 átomos,
indica que debe considerarse como el (mejor) valor predicho por nuestra técnica. Observamos además un muy
buen acuerdo con valores de otros autores20 utilizando celdas más grandes y el código VASP24. Ambas
estimaciones, sin embargo, difieren significativamente de los (abundantes) resultados experimentales.21 Este hecho
sugiere que los cálculos de
A
Q
aquí presentados difícilmente puedan ser cuantitativamente comparados contra
experimentos. Las sutilezas que rodean los aspectos cuantitativos de este tipo de cálculos para el H, ya fueron
observadas por nosotros25 en el sistema Fe-Cr, de similar estructura cristalina que el presente.
TABLA 2: Energías de activación. (*) indica extrapolación elástica.
Por otra parte, el cálculo de
A
Q
para Zr BCC también presenta problemas. Esta fase es elásticamente inestable a 0
K,
11 12
( ) 0CC
. Ambos códigos, SIESTA y VASP (resultados propios), arrojan sin embargo valores positivos
muy pequeños de este módulo. Entendemos es por esta razón que, al introducir el H en el sitio T de celdas mayores
que 16 átomos, las mismas adquieren grandes deformaciones espurias. Resulta claro además, que las
extrapolaciones elásticas en estas circunstancias serían de validez dudosa.
Resumiendo, en virtud de los dos párrafos previos decidimos, para todos los cálculos de barreras de energía que
siguen, fijar el tamaño de celda en 16 átomos y no aplicar correcciones elásticas. De este modo trabajamos con un
conjunto de valores comparables entre sí, con el cual pretendemos al menos reproducir correctamente
comportamientos cualitativos.
Aleaciones ordenadas: barreras de energía
Geométricamente podemos identificar las aleaciones estudiadas,
/16n
, mediante las coordenadas de la especie
minoritaria dentro de la supercelda de
222
. En unidades de parámetro de red,
1/16 y 15/16: (0,0,0)
2/16 y 14/16: (0,0,0) (1,1,1)
4/16 y 12/16: (0.5,0.5,0.5) (1.5,1.5,0.5) (1.5,0.5,1.5)
(0.5,1.5,1.5)
8/16: como en estructura B2.
A fin de revelar posibles sistemáticas y determinar la posibilidad de su aplicación a aleaciones aleatorias,
caracterizamos las barreras atendiendo a los octaedros (entorno cercano) centrados en sitios O, dentro de los cuales
ocurren los saltos. Para las aleaciones consideradas existen 5 de estos octaedros, que comprenden un total de 8
saltos, identificados como
n
J
y representados esquemáticamente en la Fig. 8.
FIG. 8: Los 8 saltos posibles en aleaciones ordenadas de acuerdo al octaedro continente.
Las figuras muestran una proyección en planta del sitio O de la red BCC. Los círculos pequeños indican sitios T,
sobre una cara de la celdilla unitaria corriente. Los círculos más grandes representan átomos de Nb o Zr, aquellos
en el centro quedan por arriba y por debajo del plano del papel; en (b) esos 2 átomos son distintos y en (d) iguales.
Finalmente, la especie minoritaria aparece en gris.
Las Tablas 3 y 4 reúnen todos nuestros cálculos de barreras. La barrera asociada a un dado salto
n
J
no es
necesariamente simétrica respecto del sentido del salto. A los fines de estas tablas tomamos el valor medio de
ambos sentidos. Para el cálculo del coeficiente de difusión en cambio (e.g. via KMC), tales asimetrías son
consideradas. En todos los casos en que un dado
n
J
se repite, los saltos respectivos son equivalentes desde el punto
de vista del octaedro, siempre y cuando no se cruce la composición del 50%; en este último caso los octaedros
asociados al salto
6
J
resultan complementarios (Fig. 8(d)).
TABLA 3: Barreras de energía
ΔE
(eV) según la concentración de Nb en la aleación.
TABLA 4: Barreras de energía
ΔE
(eV) según la concentración de Nb en la aleación.
De la inspección de las tablas surge claro un comportamiento sistemático de las energías de salto, como así
también que el entorno cercano definido por el octaedro resulta algo deficiente para la predicción cuantitativa. En
general las aleaciones ricas en Nb muestran un comportamiento más consistente a lo largo de las distintas
concentraciones, siendo más notable la dispersión de valores para las ricas en Zr. Nótese en particular lo que
sucede con
1
J
, donde las
E
para las aleaciones ricas en Nb prácticamente coinciden con el valor de Nb puro,
Tabla 2, no siendo así para las ricas en Zr. Queda para la especulación si este resultado está o no relacionado con la
(casi) inestabilidad elástica del Zr BCC mencionada anteriormente.
Aleaciones ordenadas: coeficientes de difusión
Con las barreras disponibles nos abocamos al cálculo de los coeficientes de difusión,
H
D
, en particular
determinamos una energía de activación efectiva,
A
Q
, por ajuste de Arrhenius en el intervalo entre 300 K y 1000
K. Los resultados de este procedimiento se consignan en la Tabla 5. Para algunas aleaciones el cálculo puede
hacerse analíticamente, en otros casos utilizamos KMC, condición que también se indica en la tabla. A modo de
ejemplo, a continuación detallamos un caso de KMC y otro analítico de cierto interés.
TABLA 5: Energías de activación (eV) entre 300 y 1000 K para las aleaciones estudiadas.
Para la composición 15/16 se identifican 4 sitios T no equivalentes,
, , ,A B C D
, en la Fig. 9, de acuerdo a su
posición respecto de un átomo de Zr. La figura también indica todos los saltos no equivalentes, en el sentido en que
ingresan al código KMC. Nótese que sitios equivalentes pueden estar conectados por saltos no equivalentes. Por
otra parte, la frecuencia del salto inverso queda determinada por la condición de balance detallado, dado las
energías relativas de los sitios T involucrados. Para este caso particular, los parámetros más relevantes del cálculo
fueron,
• Frecuencia de ataque:
13
02.396 10
Hz
• Matriz de barreras (eV):
1 2 1 3
1 2 1 3
00 0.1774 0.1022 0 0 0.2524
0 0 0 0 0.1726 0.1963 0 0
0 0 0 0 0 0.1804 0.2430 0
0 0 0 0 0 0 0.2077 0.1776
A A AB A A
BB BC
CC CD
D D D D
(6)
• Número de trayectorias: 10000
• Número de saltos/trayectoria: 10000
Los 2 últimos aseguran un error en
H
D
inferior a 1%.
FIG. 9: Saltos para la composición 15/16.
La composición 4/16 admite resolución analítica. La estructura cristalina consiste de celdillas BCC unitarias
conteniendo Nb en su centro, rodeadas en sus 6 caras por celdillas conteniendo Zr (mayoritario) y viceversa. Esa
simetría produce que los sitios T sean todos equivalentes, y los saltos entre ellos simétricos:
27
,JJ
y
8
J
de Tabla
4. De este modo identificamos caminos de migración independientes, paralelos a los ejes Cartesianos, que pueden
interpretarse como una sucesión de un paralelo de saltos
2
J
en serie con un paralelo entre
7
J
y
8
J
, representados
en la Fig. 10. Luego, aplicando ideas similares a las del modelo fenomenológico, obtenemos una expresión cerrada
para
H
D
,
2
0 2 7 8
1 12 1 1 ,
2
H
DA
(7)
donde
0
A
es el parámetro de red de la aleación y
i
las frecuencias de salto. Incidentalmente, para esta
composición
2
J
resulta un cortocircuito relativo, y siendo la barrera de
7
J
muy alta,
H
D
queda dominado por
8
J
.
De ahí el pozo en
A
Q
respecto de las concentraciones vecinas, c.f. Tabla 5. Obsérvese que para la aleación espejo
en cambio domina
2
J
.
FIG. 10: Caminos de migración en la aleación 4/16.
Los gráficos de Arrhenius
H
logD
vs.
1/T
presentan una leve desviación respecto del comportamiento lineal, no
obstante es posible identificar los saltos cuyas barreras dominan
A
Q
. En la Fig. 11 reunimos todas las barreras y
energías de activación calculadas. Aparte de los casos
2
J
y
8
J
recién mencionados, observamos que
1
J
es
dominante para las composiciones de Nb y Zr puros (obviamente), 14/16 y 15/16;
4
J
domina en las aleaciones
1/16 y 2/16; finalmente
6
J
domina para 50%.
FIG. 11: Energías de salto y de activación en las aleaciones estudiadas.
Por último ya mencionamos que no esperamos acuerdo cuantitativo de nuestros cálculos con los escasos valores
experimentales, e.g.
0.26
A
Q
eV para Zr-
20%Nb
,4 o
0.28
A
Q
eV para la fase bifásica de 31% Zr-
enriquecido +69% Zr-
.1 No obstante, en la Fig. 12 mostramos nuestros resultados para
H
D
vs. % de Nb para
diversas temperaturas. Este gráfico resulta similar a otro calculado por Samin et al.26 en el sistema Fe-Cr, BCC
desordenado, utilizando otras metodologías (mucho más costosas). En literatura se ha notado1 que a medida que la
fase
se enriquece en Nb, disminuye el volumen de la misma y eventualmente pierde continuidad, lo cual
produce un descenso de
H
D
. La Fig. 12 sugiere que dicho descenso también ocurriría de modo intrínseco, por lo
menos hasta alcanzar 50% de Nb.
FIG. 12: Coeficiente de difusión
H
D
en función del contenido de Nb y de la temperatura.
IV. CONCLUSIONES
Brevemente, desarrollamos un modelo de difusión fenomenológico adaptado de otros existentes en la literatura y
calculamos el coeficiente de difusión del H en aleaciones de Zr-Nb para diferentes concentraciones de Nb.
Respecto de la primera pregunta que planteáramos al inicio en I, podemos decir que los resultados de nuestro
modelo de difusión son consistentes con la hipótesis de que la disminución observada en la difusividad del H a
medida que transcurre el tiempo se debe a la pérdida de continuidad de las láminas de fase
que originalmente
rodean a los granos
del tubo de presión.
El modelo también arroja valores de la relación axial/radial para la velocidad de propagación de fisuras,
comparables a valores de literatura, dentro de un marco teórico mayormente aceptado, que la relacionan con la
difusividad. Esto sería una consecuencia de la textura del tubo heredada del proceso de fabricación.
Respecto de la segunda pregunta, nuestros cálculos ab initio predicen que efectivamente la difusividad del H se ve
afectada por el contenido de Nb en la fase
, y que varía de manera no monótona pasando por un mínimo para la
concentración de 50%Nb. Esto podría convertirse en una razón adicional de la disminución de la difusividad antes
mencionada. Identificamos aquí dos caminos posibles para profundizar el estudio de la cuestión: i) extender los
cálculos a aleaciones realmente desordenadas y ii) estimar el efecto de los bordes de grano en términos de tamaño
y temperatura a la manera de.27
AGRADECIMIENTOS
Queremos agradecer a la Dra. Gladys Domizzi por sus comentarios sobre el manuscrito. Este trabajo fue
parcialmente financiado por el PIP CONICET: 11220170100021CO.
REFERENCIAS
[1] B. Skinner y R. Dutton. en Hydrogen effects on material behavior (eds. Moody, R. y Thompson,W.) 73-83 (The
Minerals, Metals & Materials Society, 1990).
[2] D. Khatamian, A. Shaddick y V. Urbanic. Influence of neutron irradiation on H diffusion in Zr–2.5Nb alloy. J.
Alloys Compd. 293-295, 324-328 (1999).
[3] C. Garcia. Solubilidad de hidrogeno en aleaciones base Zr-Nb para aplicaciones nucleares Tesis de mtria.
(UNSAM, 2019).
[4] A. Sawatzky, G. Ledoux, R. Tough y C. Cann. en Metal–Hydrogen Systems (ed. Vizeroglu, T.) 109-120
(1982).
[5] J. Mieza, G. Vigna y G. Domizzi. Evaluation of variablesaffecting crack propagation by Delayed Hydride
Cracking in Zr–2.5Nb with different heat treatments. J. Nucl. Mater. 411, 150-159 (2011).
[6] G. H. Vineyard. Frequency factors and isotope effects in solid state rate processes. J. Phys. Chem. Solids 3,
121-127 (1957).
[7] J. M. Soler, E. Artacho, J. D. Gale, A. Garcia, J. Junquera, P. Ordejon y D. Sanchez-Portal. The SIESTA
method forab initioorder-Nmaterials simulation. J. Phys. Condens. Matter 14, 2745-2779 (2002).
[8] J. Kearns. Diffusion coefficient of hydrogen in alpha zirconium, Zircaloy-2 and Zircaloy-4. J. Nucl. Mater. 43,
330-338 (1972).
[9] Effects of Radiation on Materials: Sixteenth International Symposium (eds. Kumar, A., Gelles, D., Nanstad, R.
y Little, T.) (ASTM International, 1994).
[10] R. Dutton, K. Nuttall, M. P. Puls y L. A. Simpson. Mechanisms of hydrogen induced delayed cracking in
hydride forming materials. Metall. Mater. Trans. A 8, 1553-1562 (1977).
[11] S. Hanlon. The Effect of Testing Direction on DHC Growth Rate Using Zr-2.5Nb Plate Tesis doct. (Carleton
University, Canada, 2013).
[12] W. M. Young y E. W. Elcock. Monte Carlo studies of vacancy migration in binary ordered alloys: I. Proc.
Phys. Soc. 89, 735-746 (1966).
[13] G. Simonelli. Simulacion por computadora de defectos puntuales en Fe, Mo y Cr Tesis doct. (UBA,
Argentina, 1998).
[14] M. Jovanovi´c, R. Eadie, Y. Ma, M. Anderson, S. Sagat y V. Perovi´c. The effect of annealing on hardness,
microstructure and delayed hydride cracking in Zr–2.5Nb pressure tube material. Mater. Charact. 47, 259-268
(2001).
[15] Canadian Standards Association. CSA-N285.15, Technical requirements for in-service evaluation of
zirconium alloys pressure tubes in CANDU reactors (2015).
[16] Kharchenko y Kharchenko. Ab-initio calculations for the structural properties of Zr-Nb alloys. Condens.
Matter Phys. 16, 13801 (2013).
[17] P. Blaha, K. Schwarz, P. Sorantin y S. Trickey. Fullpotential, linearized augmented plane wave programs for
crystalline systems. Comput. Phys. Commun. 59, 399-415 (1990).
[18] G. Grad, J. Pieres, A. F. Guillermet, G. Cuello, J. Granada y R. Mayer. Systematics of lattice parameters and
bonding distances of the omega phase in ZrNb alloys. Phys. B: Condens. Matter 213-214, 433-435 (1995).
[19] G. Benites, A. F. Guillermet, G. Cuello y J. Campo. Structural properties of metastable phases in Zr–Nb
alloys. J. Alloys Compd. 299, 183-188 (2000).
[20] W. Lu, A. Gao, Y. Liu y Z. Dai. Diffusion behaviors of hydrogen isotopes in niobium from first-principles.
Sci. China: Phys. Mech. Astron. 55, 2378-2382 (2012).
[21] G. Schaumann, J. Volki y G. Alefeld. The diffusion coefficients of hydrogen and deuterium in vanadium,
niobium, and tantalum by gorsky-effect measurements. physica status solidi (b) 42, 401-413 (1970).
[22] C. Varvenne, F. Bruneval, M.-C. Marinica y E. Clouet. Point defect modeling in materials: Couplingab initio
and elasticity approaches. Phys. Rev. B 88 (2013).
[23] R. C. Pasianot. On the determination of defect dipoles from atomistic simulations using periodic boundary
conditions. Philos. Mag. Lett. 96, 447-453 (2016).
[24] G. Kresse y J. Furthmuller. Efficient iterative schemes forab initiototal-energy calculations using a plane-wave
basis set. Phys. Rev. B 54, 11169-11186 (1996).
[25] P. Bruzzoni y R. Pasianot. A DFT study of H solubility and diffusion in the Fe-Cr system. Comput. Mater.
Sci. 154, 243-250 (2018).
[26] A. J. Samin, D. A. Andersson, E. F. Holby y B. P. Uberuaga. First-principles localized cluster expansion study
of the kinetics of hydrogen diffusion in homogeneous and heterogeneous Fe-Cr alloys. Phys. Rev. B 99, 014110
(2019).
[27] V. P. Ramunni, M. I. Pascuet, N. Castin y A. M. Rivas. The influence of grain size on the hydrogen diffusion
in bcc Fe. Comput. Mater. Sci. 188, 110146-110156 (2021).
