Theoretical modelling of iron nitriding coupled with a nanocrystallisation treatment. Application to numerical predictions for ferritic stainless steels

Applied Surface Science 258 (2012) 6611–6620

Contents lists available at SciVerse ScienceDirect

Applied Surface Science journal homepage: www.elsevier.com/locate/apsusc

Theoretical modelling of iron nitriding coupled with a nanocrystallisation treatment. Application to numerical predictions for ferritic stainless steels B. Panicaud ∗ , M. Chemkhi, A. Roos, D. Retraint ICD-LASMIS, Université de Technologie de Troyes (UTT), UMR CNRS 6279, 12 rue Marie Curie, 10010 Troyes, France

a r t i c l e

i n f o

Article history: Received 31 January 2012 Received in revised form 14 March 2012 Accepted 15 March 2012 Available online 23 March 2012 Keywords: Nitriding Scale transition methods Nanocrystallisation Diffusion Stress effects

a b s t r a c t This paper analyses a recently developed duplex process combining nitriding with nanocrystallisation. A model is proposed to show how nitrogen diffusion mechanisms are modified within ferritic steels due to the nanostructure near the top surface. This model is based on micro-mechanical and micro-physical approaches, and also on the thermodynamics of irreversible processes. It takes into account size effects influencing the nitrogen diffusion, including mechanical stresses at the different length scales. Several models are investigated and numerical applications are performed. The results are compared to literature in order to demonstrate the generality of the present methodology. © 2012 Elsevier B.V. All rights reserved.

1. Introduction The surface mechanical attrition treatment (SMAT) is a recently developed process for nanocrystallising the surface of metallic alloys, thereby improving several mechanical properties. This has already been demonstrated for metallic materials, for instance for alloys such as steels [1,2] or pure metals such as copper [3,4]. Similarly, nitriding can improve hardness near the top surface of a material [5]. It is now well established that this is caused by diffusion of nitrogen from the surface into the bulk material [6]. However, when combining both treatments into a duplex process, the interaction between them is not really clear [7,8]. The present paper presents a model for the nitrogen diffusion in such a duplex process. Diffusion paths are expected to be modified by a smaller grain size. Furthermore, there should also be a size effect with respect to the grain size, as will be discussed in Section 2.6. Also, a trapping effect is often used to explain modification of the diffusion occurring for microstructures [9]. However, this trapping effect is not considered here, because the present diffusion model of nitrogen does not need neither sinks nor sources. The structure of the present paper is as follows: the first part presents the model of the diffusional process and its general foundations. The model is derived from the Fick equations in a multi-layer system. Its solution provides the chemical phases present in the material, as well as the extent of the layers in which

∗ Corresponding author. Tel.: +33 3 25 71 80 61; fax: +33 3 25 71 56 75. E-mail address: [email protected] (B. Panicaud). 0169-4332/$ – see front matter © 2012 Elsevier B.V. All rights reserved. doi:10.1016/j.apsusc.2012.03.089

they develop. Particular aspects of the model with and without SMAT are discussed. Then, a discretised version of the model is set up. The main hypotheses as well as the initial and boundary conditions are given. The solution of this numerical model can be obtained in different manners, each requiring different simplifying hypotheses, and the resulting numerical algorithm is described. The second part gives the numerical results, and the effects of SMAT on the nitrogen diffusion are demonstrated. Finally, some results are compared to experimental ones from the literature. Throughout the paper, the model is applied to ferritic stainless steels that have been described in the literature (for instance in [10]). By comparing experimental results and modelling, it was already shown [11] that the improvement of the mechanical properties of SMATed materials is retained throughout the nitriding process. However, those models took into account the grain size variation only through simple relations [11]. The objective here is to use scale transition methods in order to propose an improved model which includes mechanical effects. 2. Theoretical model This section presents the different parts of the model: first the chemical phases are presented from a thermodynamic point of view, as well as the general principles and assumptions that drive the model, leading to the diffusion equations and their conditions. Then, a subsequent subsection deals with the SMAT effect, as well as the numerical algorithm for solving the model equations.

6612

B. Panicaud et al. / Applied Surface Science 258 (2012) 6611–6620

Nomenclature X − → X X X X·Y X:Y X⊗Y T

scalar variable rank one tensor = vector second rank tensor fourth rank tensor contraction between tensors X and Y double contraction between tensors X and Y tensorial product between tensors X and Y ¯

T

¯ or (X) (X) transpose of X (second or fourth rank tensor)   ¯ X¯  Euclidean norm of a second rank tensor X¯

2.1. Phase diagram The binary Fe–N diagram (Fig. 1) contains three particular chemical phases ␣-Fe(N), ␥ -Fe4 N and ␧-Fe2–3 N that develop during the nitriding process. ␣-Fe(N) is a solid solution of nitrogen in ferritic steel. The phases have an extended composition domain (expressed here in weight fraction (wt%); hereafter nitrogen content corresponds to normalised weight fractions. The nitriding is carried out at a maximum temperature of 425 ◦ C, similar to what is applied usually to steels [5]. This maximum temperature sets a reasonable process duration while avoiding re-crystallisation (the latter would undo the benefits of the SMAT process). Depending on the duration of the nitriding process, the thicknesses of the layers formed by

Fig. 1. Nitrogen-iron binary phase diagram.

Fig. 2. Geometry of the ferritic steel substrate with it diffusion layers.

these phases vary between several to several hundreds of microns. One main objective of the present work is to determine these thicknesses. Three treatment durations will be modelled (5, 10 and 20 h) as well as three temperatures (350, 400 and 425 ◦ C). 2.2. General principles and hypotheses The growth of iron nitrides takes place according to the geometry shown in Fig. 2: the problem is a three-layer problem. Inside each layer, the diffusion process has to be treated at a microscopic length scale in order to take into account accurately the preferential diffusion paths such as grain boundaries. This also allows taking into account the influence of the stress at each length scale. The present model therefore combines thermodynamics of irreversible and reversible processes with scale transition methods. Its simplifying hypotheses are the following: – three length scales are distinguished: the microscopic scale, consisting of either grain (lattice) or grain boundaries (each one indicated by a superscript ); the mesoscopic scale, consisting of aggregates of grains in each of the three layers (where each layer, indicated by a subscript “i”, consists of one single chemical phase); and the macroscopic scale, which is the scale of the entire system. The latter scale is only used for the computing step, so the following sections will deal almost exclusively with the scale transitions between the mesoscopic and microscopic length scales. Thermodynamic quantities at the mesoscopic scale, such as concentration, chemical potential, stress and strain, will be indicated with an upper case notation, whereas the corresponding ones at the microscopic scale will be indicated by a lower case notation. The quantities used at the different length scales are summarised in Table 1. – volume forces (such as gravitation etc.) have no significant influence here on nitrogen diffusion; – neither sources nor sinks are incorporated in the layers; this is justified if no phase transformations occur away from the interfaces during the diffusion process;

Table 1 Summary of the different quantities depending on the considered scale. Diffusional (concentration, flux, chemical potential, other) Mesoscopic

Microscopic

c˜i˝ , c˜ilat , c˜i − →˝ − →lat − →gb j i , j i , j i

˙, ˙ i

gb

gb

˝ , lat , i i i

˝ lat gb D0,i , D0,i , D0,i gb lat ˝ Q0,i , Q0,i , Q0,i

Without scale

Mechanical (stress, pressure, volume variation) nit

˜ C˜ i , C˜ xt C, − → Jx , Ji Mi Deff,i , ˛tx

˝

lat

Geometrical

Other

hi

gb

gb

i , i , i

fi˝ , filat , fi

˝,nit lat,nit gb,nit i , i , i ˝,dif lat,dif gb,dif i , i , i gb lat ˝ pi , pi , pi ˝,res i (x, t = 0) v˝ N,i

˚lat i

T

B. Panicaud et al. / Applied Surface Science 258 (2012) 6611–6620 Table 2 Existence domains of the Fex Ny phases, in normalised concentration.

– nitrides are supposed to growth perpendicularly to the initial surface of specimens, so that symmetry of the system leads to a one-dimensional problem; – only growth of iron nitrides are considered; precipitation of other phases (such as CrN) or interaction of N with other chemical elements are not taken into account; – at the mesoscale, the linearly elastic mechanical response of each layer “i” is supposed to be isotropic, and characterised by Young’s modulus Ei and Poisson ratio vi ; – the nitriding process induces residual stresses at the mesoscale that can roughly be calculated by using a purely elastic three-layer mechanical model (with a geometry shown in Fig. 2), considering an in-plane stress state, according to [12]:

⎛ ⎛ ⎝− ⎝ ˙init (t) = Ei∗

3 

⎞

⎞

Ej∗ hj (t)⎠ ı1i + E1∗ h1 (t)⎠

V

ε −V 



V 

⎛ ⎛ + ⎝− ⎝

Ej∗ hj (t)

where the effective Young modulus E* is E/(1 − ) and hi (t) is the thickness of each layer “i”. The latter may vary with time, thereby inducing an evolution of residual stresses during the nitriding process. ıij is the Krönecker symbol. The volume variation is due to the phase transformation and calculated simply by considering the difference of molar volume Vi with respect to the bulk (layer 3). The correspondence between layer and subscript “i” is given in Fig. 2. It can be observed that mesostresses are supposed to be uniform in each layer. 2.3. Diffusion equations

˜ ˜ T, ˙) × ∂C Deff (C, ∂x

Lower nitrogen concentration (wt%)

1 = ␧-Fe2–3 N 2 = ␥ -Fe4 N 3 = ␣-Fe(N)

– 5.75 0.08

7.35 5.50 0

crit C˜ i/i+1 =

up C˜ idown + C˜ i+1

(3)

2

⎞

⎞

Ej∗ hj (t)⎠ ı2i + E2∗ h2 (t)⎠

V

  −V˛



V˛

(1) The resulting values are reported in Tables 2 and 3. According to this hypothesis, the concentration difference in one layer is higher than the real concentration difference. If the first Fick equation would have been used directly, a thicker layer would have been obtained than in reality (for constant imposed flux and constant diffusivity). The present model therefore provides an upper limit. Nevertheless, and according to the phase diagram (Fig. 1 and Tables 2 and 3), the mean absolute error between the concentration limits used here and the real ones is quite small. 2.4. Problem definition (boundary and initial conditions, grid) Initially, the system consists of one layer: the ferritic steel substrate (layer 3 in Fig. 2) and no nitrogen is present there:

Due to the nitriding process, the nitrogen concentration varies with depth (with respect to the surface), and the different Fex Ny layers can be formed; both elements constitute a Fe–N interdiffusion couple. This paragraph sets up a diffusion model of nitrogen through the ferritic steel substrate. It supposes that only nitrogen diffuses through the system. This is justified by the experimental observation that in iron nitrides Fe diffuses slower than N, as applied by several authors [10–11,13–15] for nitriding modelling. A single one-dimensional diffusion equation for a nonstationary state is considered here:



Upper nitrogen concentration (wt%)

j=1 3  j=1

∂C˜ ∂ = ∂t ∂x

Phase

with respect to the interface according to the geometry in Fig. 2), as given by the following equation:

3 

j=1

6613



(2)

where C˜ is the nitrogen concentration at mesoscopic scale, t the time and x the depth with respect to the surface; Deff is the effective mesoscopic diffusion coefficient, which varies for the different Fex Ny layers that can grow. They depend on temperature T and mesoscopic stress ˙, and Eq. (2) is therefore non-linear. Here the problem has been simplified: Deff is constant in every layer, but dif

˜ T, ˙ . However, fers between layers so that Deff = Deff layer(C), it does not take into account the concentration discontinuities at interfaces because, according to the phase diagram and the Gibbs rule, there are concentration discontinuities between each layer at mesoscopic scale. The locations of these discontinuities are determined here using comparison criteria on the interface concentrations, related to the phase diagram. This determination consists in comparing critical concentrations at the interfaces, thereby ensuring concentration continuity between each layer. These critical concentrations have been calculated as the average between the concentration limits of the lower layer “i” and the concentration limits of the upper layer “i + 1 (here upper and lower are taken

˜ C(x, t = 0) = 0,

∀x ∈ ˛ − Fe

(4)

Moreover, the two opposing faces of the sample are supposed to be identical throughout the nitriding process, so the problem is symmetrical (in reality, depends on the sample position in the nitriding apparatus). Thus, the flux must necessarily be zero in the middle of the sample: Jx (x = L, t) = 0, ∀t

(5)

where L corresponds to the half-width of the sample. The second boundary condition, i.e. the nitrogen flux at the surfaces, is more difficult to quantify because it depends on the specific nitriding process. It can be calculated from the nitrogen gas flux within the enclosure. However, approximations are necessary and the accuracy of such calculations remains limited due to the complexity of associated processes such as adsorption and dissociation. Consequently, the boundary condition at the gas/metal interface is taken as a limit concentration C˜ s of nitrogen for the substrate material, and this may also depend on time: ˜ = 0, t) = C˜ s , ∀t C(x

(6)

The formation of the nitride coatings is considered as a onedimensional semi-infinite diffusion problem. This requires that the Table 3 Critical nitrogen concentrations for the interphase positions. Interphase

Average nitrogen concentration (wt%)

␧-Fe2–3 N/␥ -Fe4 N ␥ -Fe4 N/␣-Fe(N)

6.55 2.79

6614

B. Panicaud et al. / Applied Surface Science 258 (2012) 6611–6620

diffusion depth is much smaller than the initial thickness of ferritic steel. Indeed, with this hypothesis the depth of diffusion can be roughly estimated from the relation: Ldiffusion ≈



Deff × t

(7)

For t = 20 h and Deff ≈ Dmean Ldiffusion = 85 ␮m L = 3 mm.

= 10−13

m2 /s

[16], this leads to

2.5. Model without SMAT The objective of the present paragraph is to justify Eq. (2) and to give the diffusion coefficient in each layer as an explicit function of stress and temperature. The nitrogen flux through each aggregate of the mesoscopic layer “i” is calculated using: − → Ji =



− →˝ fi˝ ji

The chemical potentials of nitrogen within each mesoscopic layer are supposed to be given by a standard relation [17]:

Mi0

Mi =

mesoscopic scale, the localisation stress tensor B is isotropic, so that its trace can be calculated as follows:

(9)

˝

where Li is the second rank Onsager tensor (the self-coupling effect of the nitrogen chemical potential with nitrogen diffusion flux), and v˝ a thermodynamic coefficient coupling a Nitrogen feature with N,i pressure p˝ . The microscopic hydrostatic pressure p˝ is calculated i i

from the trace (Tr) of the microscopic stress tensor ¯ i˝ in  of layer “i” (defined later on):



(10)

As per hypothesis, volume forces are not taken into account. Moreover, this specific dependence on stress is not the most general, but it is often used in the literature [18,19]. Chemical potential and stress have to be related to the mesoscopic solicitation calculated by a scale transition method. These methods been extensively discussed in the literature for thermodynamics and mechanics [20]. In the present application, the microscopic stress ¯ i˝ can be separated into two contributions: the effect of phase transformation due to the nitriding process and the self-induced stress effect within the microstructure due to diffusion [19]: ˝ i

=

˝,dif + i (˝ ) i

˝,nit i

˝,nit

Moreover,  i part

nit ˙ i in ˝

(11)

can be directly related to its mesoscopic counter-

each layer “i”, through a fourth rank stress localisation

tensor Bi [20]: ˝,nit

i

˝

nit

= Bi : ˙ i

(12)

Similarly, the microscopic chemical potential gradient in  of layer “i” is also related to the mesoscopic diffusion force (thermodynamically speaking), using a second rank localisation tensor

˝ Ai

[21]:

˝ − − → ˝ → ∇ i = Ai · ∇ Mi

(13)

where Mi is the mesoscopic chemical potential of nitrogen in layer “i”. Incorporating Eqs. (10)–(13) in Eq. (9) yields:



˝ − →˝ ji = −Li ·



=

B

˝ − − →˝ → − → ˝ ji = −Li · ∇ ˝ + v˝ ∇ pi i N,i

˝ Ai

− →

· ∇ Mi −

v˝ → N,i − 3

∇ Tr



˝

Bi :

nit ˙i



˝,dif + i (˝ ) i

(14)

(15)

with Mi0 and C˜ i0 respectively the chemical potential and concentration at reference pressure P0 and temperature T0 of the considered solid solution, in layer “i” at the mesoscopic scale. R is the gas constant. Because of the assumption of mechanical isotropy at the

(8)

with  corresponding either to lattice or grain boundary, fi˝ the − →˝ corresponding relative volume fraction and ji the microscopic flux in  of each aggregate in layer “i”. Each microscopic flux in  of layer “i” can be related to the corresponding microscopic chemical potential ˝ and microscopic hydrostatic stress p˝ [17]: i i



i

RT − →˜ − → ⇒ ∇ Mi = ∇ Ci C˜ i

˝

p˝ ≡ − 31 Tr ¯ i˝ i

+ RT ln

C˜ i C˜ 0

⇒



B1 ı˛ˇ ı  + B2 ı˛ ıˇ + ı˛ ıˇ





e˛ ⊗ eˇ ⊗ e ⊗ e

B : ˙

=

2B1 ˙ + B2 (Tr˙) 1

B:˙

=

(2B1 + 3B2 )Tr˙ ≡ B0 Tr˙

⇒ Tr

(16)

with 1 the second rank identity tensor. Eq. (16) can be applied indifferently in each  of layer “i”, but its numerical value is different because of different stresses in each layer. Eqs. (15) and (16) are then used to simplify Eq. (14):

v˝ nit RT ˝ ˝ − − →˝ → → N,i ˝ ˝ − ji = − Li · Ai · ∇ C˜ i + B0,i Li · ∇ Tr˙ i + ˜ 3 Ci





˝,dif (˝ ) ∂ Tr i v˝ i N,i

3

⎛ ⎝



∂˝ i

˝ − → Li · ∇ ˝ = i



˝,dif ∂ Tr i (˝ ) v˝ i N,i

3

v˝ N,i 3

∂˝ i

⎞ − 1⎠

RT ˝ ˝ − → Li · Ai · ∇ C˜ i + C˜ i

nit → ˝ ˝ − Li · ∇ Tr˙ i B0,i

(17)

The last terms of the second relation have been obtained by derivation of the diffusion stress while considering Eq. (13). The Onsager tensor L¯ can be replaced by the (microscopic) diffusion coefficient ¯ which is defined as [17]: tensor D, D≡L

RT c˜

(18)

Moreover, the diffusion tensor (or the Onsager tensor) itself also depends on stress. The coupling between both tensors is presently done at microscopic scale. For example, several authors suggested [17]:

 f + pv N N

D = D0 exp −

1

RT  f  nit vN N dif ∼ (B0 Tr˙ + Tr ) 1 1+ = D0 exp −  RT  3RT nit f v = D0 exp −

N

RT

1+

N

3RT

(B0 Tr˙

+ Tr

dif

(19)

)

fN ≡ Q0 is the Helmholtz free energy identified to an effective activation energy Q0 , and vN is the volume change in each layer occurring at the atomic scale due to the mechanical effect during the nitriding process. Note that the dependence with such an Arrhenius-like law (based on statistical Boltzmann distribution) could be also considered directly at the mesoscopic scale, depending on which input data are available.

B. Panicaud et al. / Applied Surface Science 258 (2012) 6611–6620 nit

caused by an equivalent volume variation (for a defect concentration c˜ ) can be calculated by [19]:

Recalling that ˙ i is uniform in each layer “i”, its spatial gradient vanishes in the volume of each layer and the final relation for flux diffusion can be evaluated using: − → Ji =

⎛



= −⎝

⎛



= −⎝



⎛ fi˝ ⎝1 −

v˝ N,i

∂

˝

∂˝ i

˝

⎛ ∼ = −⎝

fi˝ ⎝1 −

˝,dif Tr i

3



⎛ ⎞ ⎠

c˜i˝



C˜ i



lat,dif

− →˝ fi˝ ji

⎞

˝

3

∂˝ i



˝ D0,i

˝ · Ai

exp

−

˝ Q0,i



˝

˝

fi˝ D0,i · Ai exp

−

˝ Q0,i

RT

⎠

nit

C˜ i

1+

RT ˝,dif

lat Ei vN,i (∞) 1 + vi vlat i

(23)

⎞ c˜i˝



≡ −Deff,i layeri (C˜ i ), C˜ i , c˜i˝ , T, ˙ i ,  i



= 2˜cilat

Tr i

˝,dif (˝ ) ∂ Tr i v˝ i N,i

6615

˝

˝

D i · Ai



v˝ N,i 3RT

→˜ ⎠·− ∇ Ci

nit ˝ B0,i Tr˙ i

˝,dif

+ Tr i



⎞ → ˜ (20) ⎠·− ∇ Ci



− → (˝ ∝ c˜i˝ ) · ∇ C˜ i i

 ⎞⎞ ⎛ ˝,dif ˝ ∂ Tr  v˝ v i nit →˜ N,i ˝,dif ˝ ⎠⎠ · − ⎝1 + N,i B0,i Tr˙ i + Tr i − ∇ Ci ˝ 3RT

3

˝

The last expression is obtained by adopting some simple approximations, such as nitrogen concentration that are taken roughly of same order whatever the considered scale (∀i, c˜gb,i ∼ = c˜lat,i ∼ = C˜ i ) and developing the stress part of the expression to first order. Firstly, it can be observed that the diffusion coefficient remains a second rank tensor. This fact will be developed further. In addition, the effective behaviour of diffusion process can be observed to depend on: - the temperature T through an Arrhenius law; nit

- the mesoscopic stress ˙ due to the nitriding process induced by phase transformation; dif - the microscopic stress  due to diffusion in a self-constraint medium, linked directly to the microscopic nitrogen concentration c˜ (or the chemical potential ); - the localisation tensor A that defines the configuration of diffusion paths inside each material.

∂i

where vlat (∞) is the volume variation in each phase occurring at N,i the atomic scale due to nitrogen if it would have been embedded in an infinite lattice, and vlat the volume of the finite lattice. Eq. (23) is i ˝,dif

compatible with Eq. (11), because we have Tr i (˝ (ci˝ )). Also, i it depends linearly on the geometry of the medium. For large grains (as in untreated materials) this term vanishes, whereas for small grains (as in SMATed materials) this effect can become significant and leads to a modification of the diffusion amplitude, as for lattice: lat,dif (˚lat ) i lat,dif Tr i (˚lat ) i

noSMAT : Tr i SMAT :

→0

(24)

is the average grain size in layer “i”, roughly calculated where ˚lat i



as ≈

3

vlat . Consequently, incorporating Eqs. (21)–(24) in Eq. (20), i

the effective diffusion tensor will be modified. All effects can be studied separately through numerical simulations in order to obtain the sensitivity of each parameter.

2.6. Modelling with SMAT In this section, the model is generalised to take into account the SMAT. The latter is expected to modify the diffusion mechanisms in three ways: firstly, the spatial distribution of the grain boundaries is modified by the SMAT, as observed experimentally in steels [1,2]. This distribution influences directly the relative fractions of the lattice and of the grain boundaries: fi˝ → fi˝ (SMAT )

(21)

Secondly, the severe plastic deformation near the surface induces residual stresses [22]. These stresses operate at the microscopic scale and have to be taken into account during the nitriding process, by modifying Eq. (12) to: ˝,nit+SMAT

i

˝

nit

= Bi : ˙ i

˝,res

+ i

(x, t = 0)

(22)

Thirdly, a size effect is expected because of the nanocrystallisation treatment. This size effect can be simply understood as the deformation of a finite (as opposed to infinite) lattice or grain boundary. In the particular case of a spherical inclusion (such as a nitrolat,dif gen atom diffusing through a lattice), the associated stress  i

2.7. Simulation environment and algorithm of the diffusion model A finite-difference method has been used with a Crank–Nicholson method [23] to solve the second Fick law presented in Eq. (2). This method has the advantage of being unconditionally stable. However, the diffusion coefficient itself depends on the phase, leading to a non-linear partial differential equation (see Eq. (20)). Although commercial packages now offer effective solvers for the diffusion equations of this apparent simple form, the non-linearity of the considered diffusion coefficients demands specific attention to numerical details. Space is discretised at the mesoscale with a regular step size x according to: x =

L Nx

(25)

The most suitable number of space nodes, Nx , will be determined later from a convergence analysis. Separation between scales has also to remain valid, whatever the time. The time increment is defined as: t =

ttot Nt

(26)

where ttot represents the total nitriding process duration and Nt is the number of time increments. The duration of the nitriding treatment depends on the specific process. The influence of this parameter on the concentration profiles is presented below.

6616

B. Panicaud et al. / Applied Surface Science 258 (2012) 6611–6620

The result of the differentiation method gives an equation that t t is a function of C˜ xt , C˜ x−x , C˜ x+x (with the subscript x refers to the space node and the superscript t to the time increment), and can usually be solved a priori by a Thomas algorithm [23]. However, the diffusion coefficient must be evaluated at the midpoint t of the interval, i.e. D(C˜ x±1/2x ), which does not exist. This is why Garcia et al. [24] or Wei et al. [25] did not use implicit methods (Crank–Nicholson for instance), a priori badly suited [24,25], and they preferred to use explicit methods. Their approach imposes a constraint relation between the time increment and the space step, known as the Neumann criteria. Nevertheless, the Crank–Nicholson method is used here in order to keep the space step and the time increment independent of each other. This was already done successfully for an aluminising process [26]. In order to solve the present problem, some non-linearity needs to be taken into account to improve this Crank–Nicholson method. First of all, for each time increment, two calculations are done. The first step delivers an approximate concentration C˜ x’t+t . The second improves its precision by estimating the diffusion

t+1/2t coefficient at t + 1/2t by Deff ≈ D C˜ xt + C˜  xt+t /2 . The concentration can thus be determined again at time increment t+1/2t t . In other words, the first t + t but with Deff instead of Deff calculation determines the concentration at every time increment t using the diffusion coefficient at t − t. The second calculation yields the value of the diffusion coefficient at every space node x at the average concentration between increments t and t + t (using values at t + 1/2t) that has just been calculated. This method doubles the number of operations and increases the calculation time, but it is much more precise. Secondly, to take into account the variation in the diffusion coefficient with the space step (due to different layers and/or concentration dependence of diffusion coefficient), it is necessary to modify the basic expression of the Crank–Nicholson method. The diffusivity ˛ (defined in Eq. (27)) must be estimated at the midpoint between two successive space nodes x and x + x. Thus, this diffusivity at the half space step is interpolated linearly according to: ˛tx±1/2x

=

˛tx + ˛tx±x 2

and

˛tx

≡

t Deff (C˜ xt ) ×

t (x)

(27)

2

With this, and if the space step is sufficiently small, the nonlinear partial differential equation can be brought back into the implicit Crank–Nicholson scheme to give the following expression (Eq. (28)), which can be expressed as a “three-diagonal” matrix. The improved Crank–Nicholson method can then be used to solve this non-linear diffusion problem, using the Thomas algorithm:

˛ + ˛ x x−x ˜ t+t

−Cx−x

t+t C˜ x+x



2

+ 2C˜ xt+t

˛ + ˛ x x+x

2C˜ xt 1 −

1 2



2 ˛x +

t = C˜ x−x



1 1+ 2



+ ˛x−x ˛ ˛x + x+x 2

˛ + ˛ x x−x

˛x+x + ˛x−x 2



2 t + C˜ x+x



−

+

˛ + ˛ x x+x 2

(28)

where diffusivity ˛ is evaluated either at t, either at t + 1/2t, according to the previous remark. The unknown concentrations are on the left-hand side of this equation and the known concentrations on the right-hand side. The algorithm can be then run. First, the following quantities are initialised at t = 0: • The grain boundary and lattice distributions f ˝ ; i • The critical concentrations at interfaces C˜ crit ; i/i+1

˜ ˜ = • Initial and boundary conditions C(x, t = 0); Jx (x = L, t); C(x 0, t);

• The effective diffusion coefficients of each layer (with dependence on temperature, stresses, localisation tensor and SMAT nit

˝,dif

etc.) Deff,i T, ˙ i ,  i

˝

, Ai , SMAT ;

• The diffusivities used in the Thomas algorithm ˛tx ; Then for each time increment the values of each space node are found as follows: • For every space node:  • first calculation of the concentration C˜ xt+t ;  • calculation of an average concentration (C˜ xt + C˜ xt+t )/2; • update of the effective diffusion coefficients Dt+1/2t ; eff • calculation of the new diffusivity at this average concentration t+1/2t ˛x ; • second calculation of the true concentration with the new diffusivity C˜ xt+t . • Calculation of interface positions of every phase in the diffusion layer; • Calculation of thicknesses of every phase from the difference between the consecutive interface positions hi (t); • Update the boundary conditions (if necessary); • Calculation and/or update of the diffusion coefficients of each phase, in the diffusion layer (with dependence on temperature, stresses, localisation tensor and SMAT)  nit

˝,dif

Deff,i T, ˙ i ,  i

˝

, Ai , SMAT ;

• Calculation with update of the diffusivity used in the Thomas algorithm ˛x ; • Increment the time t → t + t. The time loop ends at ttot , the duration of the nitriding treatment. Following this algorithm, different configurations have been tested and resolved, as presented in the next section. 3. Numerical predictions 3.1. Effect of the scale transition model The first case is considered without taking into account any stresses but with explicit temperature dependence. The effective diffusion tensor reduces to: Deff,i (T ; phasei (C˜ i )) = =



 ˝

gb gb gb fi Di (T ) · Ai



˝ fi˝ D0,i

˝ · Ai

exp

lat

−

˝ Q0,i

RT



(29)

lat gb + (1 − fi )Di (T ) · Ai

˝

One more assumption is introduced. All localisation tensors Ai are supposed to be isotropic (diagonal and spherical). Considering ferritic phase (␣-Fe(N), the latter condition is directly satisfied thanks to its cubic lattice, if we suppose that distortion leading to a possible stress-induced anisotropy is a second-order approximation. The same consideration holds for ␥ -Fe4 N (cubic lattice). For ␧-Fe2–3 N (hexagonal lattice), localisation tensor depends on mesoscopic diffusion tensor that is supposed isotropic (grains are randomly distributed according hypothesis given in Section 2.2). However, it depends also on microscopic diffusion tensor linked to the crystallographic symmetry. Nevertheless, recent works have shown an influence of orientations [27], but as important as those calculated using different values for diffusion coefficients found in literature for an isotropic model. This justify that the diffusion coefficient has been considered (in first approximation) as a scalar in the previous algorithm (Section 2.7), whatever the layer “i”. Therefore, the following coefficient will be directly treated as scalar quantities, as only a one-dimensional problem is considered.

B. Panicaud et al. / Applied Surface Science 258 (2012) 6611–6620

6617

In the literature, many scale transition methods exist for obtaining the mesoscopic diffusion coefficient from the microscopic ones, such as Maxwell–Garnett, Einstein, Hart–Mortlock, etc. [28,29]. Only six of them will be used in the present paper. The first one corresponds to a parallel configuration between grain boundaries and lattice (using the arithmetic mean), called Hart’s model: ˝

// gb gb gb ∀˝, Ai = 1 ⇒ Deff,i = fi Di + (1 − fi )Dilat

(30)

Another model is the serial configuration between grain boundaries and lattice (using the harmonic mean):



˝



−1

gb gb ⊥ ∀˝, Ai → ∞ × 1 ⇒ Deff,i = fi Di



gb



+ (1 − fi ) Dilat

−1

−1

(31) In the third approach the problem is represented as an inclusion of microscopic entities with given diffusion coefficients embedded in an infinite matrix, as in the mechanical Eshelby inclusion problem [20]. This corresponds to an upper bound (similar to the Maxwell+ bound for a thermal problem or the Hashin–Shtrikman+ bound for a mechanical or magnetical equivalent problem [20]), in which the matrix is supposed to have the diffusion coefficient of the lattice: −1 (Di˝ + 2Dilat )

˝

Ai = 

fi˝ (Di˝ + 2Dilat )

−1

+ × 1 ⇒ Deff,i =

˝ gb

gb

gb

2(1 − fi )Dilat + (1 + 2fi )Di gb

gb

(32)

gb

(2 + fi ) + (1 − fi )(Di /Dilat )

The fourth estimation corresponds to a lower bound (similar to the Maxwell− bound for thermal problems or the Hashin–Shtrikman− bound for mechanical and magnetical equivalent problems [20]), in which the matrix is supposed to have the diffusion coefficient of the grain boundaries: ˝ Ai

gb −1

=

(Di˝ + 2Di )

 ˝

− ⇒ Deff,i =

gb −1

fi˝ (Di˝ + 2Di )

×1 (33)

gb gb gb 2fi Di + (3 − 2fi )Dilat gb gb gb (3 − fi ) + fi (Dilat /Di )

The fifth estimation corresponds to the self-consistent approach (SC) where the matrix is supposed to have an a priori unknown value, which can only be defined implicitly. Using the Eshelby tensor, it can be demonstrated that: ˝

Ai =

1



SC 1 + (1/3) Deff,i

−1 

SC Di˝ − Deff,i

gb



SC 1 + (1/3) Deff,i gb



−1

SC ) (Dilat − Deff,i

+

(34)

gb

fi Di

SC 1 + (1/3) Deff,i

−1

gb

SC ) (Di − Deff,i

The sixth and last model is a particular average corresponding to the geometric mean, with no particular signification for localisation tensor

˝ Ai :



× ≡ Dilat Deff,i

(1−f gb ) i



gb

× Di

f gb i

models, for two arbitrarily chosen values for the microscopic diffusion coefficients such that Dgb > Dlat , which correspond to the present conditions for most of phases (see in Table 4). In metallurgical problems, this condition would depend on temperature and SMAT treatment, and it might even be inversed, which would also invert the upper (HS+) and lower (HS−) bounds. The figure shows significant differences. The serial and parallel models respectively give the minimum and maximum values. All other models have their values in between, as would be expected. For the nitriding process with SMAT, the results of the six models are presented in Fig. 4a and a close-up in Fig. 4b. The parameters used in the calculations are given in Table 4. The values for the nitrogen penetration are clearly separated into two groups, showing the effect of grain boundaries. Contrarily, for the nitriding process without SMAT, it was observed (graphs not presented here) that the differences between the six models are very small. 3.2. Influence of the variation of the average grain size with depth Because of Eq. (21) and Eqs. (30)–(35) it can be expected that the variation of the average grain size with depth affects the results, specifically after the nanocrystallisation due to the SMAT. Literature provides different variations of the average grain size with depth for steels [15]. On one hand, the simplest one is considering the distribution as a succession of three layers: 1. a top layer with nano-grains; 2. a transition layer with a linear transition from nano-grains towards micro-grains with increasing depth; 3. a core layer with unmodified grains of micrometer size.

×1

(1 − fi )Dilat

SC = ⇒ Deff,i

Fig. 3. Effective diffusion coefficient using the six scale transition models.

(35)

As an example, in Fig. 3 the mesoscopic diffusion coefficients are plotted as a function of relative grain boundary fraction for all six

Studies on ferritic steels give their layer thicknesses, as summarised in Table 5. On the other hand, research on steels [1,2] has shown that the variation of the average effective grain size could be represented more precisely using the following relation (hereafter called logarithmic variation), based on the mean free path of dislocations: ˚(x) =



ϕ0 ϕ1 − ln x

2

, ∀x < xc

(36)

where ϕ0 and ϕ1 are parameters and xc is a maximum depth at which the average grain size becomes the one of the unmodified core layer. The parameter values for such a distribution are also given in Table 5, and both average grain size variations are shown in Fig. 5.

6618

B. Panicaud et al. / Applied Surface Science 258 (2012) 6611–6620

Table 4 Values for interdiffusion coefficients of the different Fex Ny phases found in literature [11–15]. * = values for activation energies and pre-exponential coefficients are not provided; only the effective diffusion coefficients are directly given from literature (corresponding to estimated values at 300 ◦ C). For given activation energies and preexponential coefficients, the diffusion coefficients are calculated at 425 ◦ C. Phase

Micro-structure

Pre-exponential factor D0 (m2 /s)

Activation energy Q0 (J/mol)

Diffusion coefficient D (m2 /s)

␣-Fe(N) ␣-Fe(N) ␥ -Fe4 N ␥ -Fe4 N ␧-Fe2–3 N ␧-Fe2–3 N

Lattice Grain boundary Lattice Grain boundary Lattice Grain boundary

6.60 × 10−7 * 1.68 × 10−9 * 2.10 × 10−8 *

77,900 * 64,000 * 93,517 *

9.8 × 10−13 3.67 × 10−11 2.73 × 10−14 2.95 × 10−15 2.12 × 10−15 2.87 × 10−13

Fig. 5. Variation of the average grain size versus depth for different distribution models.

from experiment is needed to have a correct prediction for nitrogen diffusion. Conversely, for a given measurement of the concentration profile this could be used as an inverse method to obtain information on the variation of the average grain size with depth.

Fig. 4. Nitrogen contents versus depth for different scale transition models without any stress dependence (linear SMAT distribution; Lgb = 0.5 nm; ttot = 20 h; T = 425 ◦ C). (a) is the global plot; (b) is a magnification part. It can be observed that the models are separated into two groups. The geometrical model and the self-consistent one are superimposed.

The effect on nitrogen concentration has also been computed and is presented in Fig. 6, which shows that the average grain size distribution after SMAT leads to an effect on predicted nitrogen content which can be small if considering a simple distribution, but which can be quite significant for a more realistic one. Indeed, based on strong hypothesis on grain size variation, the influence of its mathematic dependence is significant. This means that an accurate measurement of such a variation of the average grain size

Fig. 6. Nitrogen contents versus depth for different variations of the average grain size without any stress dependence (parallel/Hart’s model; Lgb = 0.5 nm; ttot = 20 h; T = 425 ◦ C).

Table 5 Parameters for both variation of the average grain size with depth: linear SMAT and logarithmic SMAT. Linear SMAT Grain size Top layer

10 nm

Transition layer Core layer

Linear variation from 10 nm to 25 ␮m 25 ␮m

Logarithmic SMAT Depth 0–15 ␮m

15–200 ␮m >200 ␮m

Grain size

Depth (xc )

Calculated with ϕ0 = 18.72 ϕ1 = 5.39 – 25 ␮m

0–200 ␮m

– >200 ␮m

B. Panicaud et al. / Applied Surface Science 258 (2012) 6611–6620

Fig. 7. Nitrogen contents versus depth for different grain boundary thicknesses without any stress dependence (parallel/Hart’s model; linear SMAT distribution; ttot = 20 h; T = 425 ◦ C). The curve without SMAT is plotted for only one grain boundary thickness.

Another parameter has to be taken into account: the grain boundary thickness. Each scale transition (Eqs. (30)–(35)) is computed with a lattice fraction that depends directly on grain size, but gb also on grain boundary thickness as a definition of fi . The present model supposes that the grain boundary thickness is homogeneous and independent of the process parameters or of the position within the sample. Consequently, several grain boundary thicknesses have been tested in order to show their influence. The results, computed with the parallel model of Eq. (30), are presented in Fig. 7. The effect of SMAT is calculated for different grain boundary thicknesses. For a given depth, it shows that the nitrogen content is higher for thicker grain boundaries, equivalent to a larger relative grain boundary fraction for a given grain size. The difference between the curves for 0.5 and 50 nm is only 25%. Nevertheless, it indicates that this parameter is significant and should be determined accurately to improve the precision of the numerical model. 3.4. Influence of the mesoscopic stress The dependence of the diffusion coefficient on the mesoscopic stress in Eq. (22) has also to be taken into account. The coefficient v˝ B˝ is assumed to be the same in the lattices and in the grain N,i 0,i boundaries (whatever ). The microscopic stress is neglected in this calculation. The effective diffusion tensor now becomes: nit

nit



×

1+

vN,i B0,i nit Tr˙ i 3RT nit

where Deff,i (T, ˙ i



Fig. 8. Nitrogen contents versus depth for mesoscopic stress dependence but no microscopic stress (parallel/Hart’s model; linear SMAT distribution; Lgb = 0.5 nm; ttot = 20 h; T = 425 ◦ C). Table 6 Different parameters for numerical simulations: coupling terms of diffusion parameters with mesoscopic stresses, and number of spatial and time steps used for calculations.

3.3. Influence of the grain boundary thickness

Deff,i (T, ˙ i ; phasei (C˜ i )) ∼ = Deff,i (T, ˙ i

6619

Parameter

Value

Coupling term with mesostress Linear coupling term with microstress Nx Nt

1.081 GPa−1 3.880 1000 36,000

meso-stress in diffusion coefficient of Eq. (37), which is around 1.1 GPa−1 [19], as summarised in Table 6. 3.5. Effect of the nitriding temperature and nitriding time Finally, the effect of the nitriding temperature on the diffusion process is calculated. The results are shown in Fig. 9 for 350, 400 and 425 ◦ C. In addition, the influence of nitriding time (5, 10 and 20 h) is also investigated and presented in Fig. 10. As expected, time and temperature have a similar effect: the higher the temperature or the time, the higher the nitrogen content at a given depth. However, the temperature or time have a limited influence as compared to other parameters. For instance, for a temperature variation T of 75 ◦ C or a time ttot multiplied by four, the nitrogen content only changes by 33% at a depth of 75 ␮m. Remark that higher values of the temperature variation would have given the same trend.

= 0 ; phasei (C˜ i )) (37)

= 0 ; phasei (C˜ i )) is given by Eq. (29).

nit The stress ˙ i

in each layer “i” should be calculated properly by using Eq. (1) for all time steps. However, to show the influence of this stress in a more direct manner, and to obtain a first indication of the order of its magnitude, two different mesoscopic stresses are imposed in the present study: 100 MPa and 1000 MPa. The results are presented in Fig. 8. For 1000 MPa, the influence would be significant, with a maximum increase of nitrogen content of 20% at a depth of around 100 ␮m. Numerical results show that the influence is more pronounced in phase ␥ than in the other phases, because of the “high” value of its diffusion coefficient. For 100 MPa, the effect is not significant. This is essentially due to the coupling value of

Fig. 9. Nitrogen contents versus depth for different nitriding temperatures without any stress dependence (parallel/Hart’s model; linear SMAT distribution; Lgb = 0.5 nm; ttot = 20 h).

6620

B. Panicaud et al. / Applied Surface Science 258 (2012) 6611–6620

Fig. 10. Nitrogen contents versus depth for different nitriding times without any stress dependence (parallel/Hart’s model; linear SMAT distribution; Lgb = 0.5 nm; T = 425 ◦ C).

Consequently, other parameters can be as significant as temperature or time. 3.6. Discussion on numerical results Comparison of numerical results on ferritic steels can be done with experimental results as presented in [11]. There is a good qualitative agreement between both approaches. Indeed, the shape of the calculated nitrogen content profile corresponds to the experimental ones. However, discrepancies between diffusion parameters lead to a difference between the calculated diffusion length and the experimental one. Especially for ␥ phase, literature proposes different values that can explain the important length of the corresponding layer with its chosen value in the present case. Despite this dispersion, the results on the concentration profile are correct. Thus, the precision of the model is acceptable in comparison with the uncertainties of experimental results, justifying the assumptions made. It is also important to note that some other approaches in ferritic steels has been performed (see for instance [30]) that corresponds to more or less the same methodology and assumptions, but including other phenomena such as precipitation. In the present study, details on micro-mechanical and micro-physical dependence of parameters is the originality of this work. 4. Conclusions In the present paper, a multiscale and multiphysics model is proposed to describe the nitriding process of ferritic steels. With reasonable approximations, it gives an expression for the mesoscopic diffusion tensor, and when coupled with a suitable discretisation in space and time of the diffusion equations, nitrogen concentration profiles can be calculated for different cases. The following geometrical, mechanical and thermal aspects can be taken into account: the diffusion path through the microstructure, the temperature, the mesostress and the microstress. The effect of microstress, which is mainly due to diffusion, can be significant for example in austenitic steels and remains yet to investigated. The paper also demonstrates that for the prediction of the nitrogen concentration profiles for the duplex process of SMAT followed by nitriding, the mechanical and morphological modifications

induced in the microstructure by SMAT can and should be taken into account as accurate as possible. These modifications influence the relative grain boundary fraction, the initial distribution of residual stresses and a size effect due to diffusion within small finite geometries such as nano-grains. The influence of these parameters has been tested numerically. Instead of classical Hart’s model, several other scale transition models have been used, depending on the assumed topology of the microstructure. However, the numerical results show that the results of all models roughly fall into two separate groups. For example, the self-consistent approach gives similar results as the parallel/Hart’s model. As expected, for a given model, varying the grain boundary thickness and the variation of the average grain size with depth leads to different results, because they directly influence the grain boundary fraction. However, it seems that the initial variation of average grain size with depth is of particular importance to obtain accurate predictions. The influence of mesostresses has also been calculated. Their influence can be of the same order of magnitude as the temperature or time parameters. This means that the initial metallurgical state as well as its evolution is of great importance to predict nitriding process, especially with a refined microstructure. References [1] T. Roland, D. Retraint, K. Lu, J. Lu, A445, Mater. Sci. Eng. (2007) 281–288. [2] L. Waltz, D. Retraint, A. Roos, P. Olier, J. Lu, Mater. Sci. Forum 614 (2009) 249–254. [3] R. Blonde, H.-L. Chan, N. Allain-Bonasso, B. Bolle, T. Grosdidier, J. Lu, J. Alloys Compounds 504 (2010) 410–413. [4] Y.M. Wang, K. Wang, D. Pan, K. Lu, K.J. Hemker, E. Ma, Scripta Mater. 48 (2003) 1581–1586. [5] Source Book on Nitriding, ASM, Metals Park, OH, 1977. [6] T. Liapina, A. Leineweber, E.J. Mittemeijer, Scripta Mater. 48 (2003) 1643–1648. [7] W.P. Tong, C.Z. Liu, W. Wang, N.R. Tao, Z.B. Wang, L. Zuo, J.C. He, Scripta Mater. 57 (2007) 533–536. [8] W.P. Tong, N.R. Tao, Z.B. Wang, H.W. Zhang, J. Lu, K. Lu, Scripta Mater. 50 (2004) 647–650. [9] S. Parascandola, W. Moller, D.L. Williamson, Appl. Phys. Lett. 76 (2000) 16. [10] I. Campos, R. Torres, O. Bautista, G. Ramırez, L. Zuniga, Appl. Surf. Sci. 249 (2005) 54–59. [11] M. Keddam, M.E. Djeghlal, L. Barrallier, Appl. Surf. Sci. 242 (2005) 369–374. [12] S. Timoshenko, J.N. Goodier, Théorie de l’élasticité, Librairie Polytechnique, Paris, 1961. [13] M. Keddam, M.E. Djeghlal, L. Barrallier, Mater. Sci. Eng. A 378 (2004) 475–478. [14] T. Belmonte, M. Goune, H. Michel, Mater. Sci. Eng. A 302 (2001) 246–257. [15] P. Cavaliere, G. Zavarise, M. Perillo, Comp. Mater. Sci. 46 (2009) 26–35. [16] H. Wang, W. Yang, Scripta Mater. 50 (2004) 529–532. [17] S.R. De Groot, P. Mazur, Non Equilibrium Thermodynamics, Dover Publications, New York, 1984. [18] A. Galdikas, T. Moskalioviene, Comp. Mater. Sci. 50 (2010) 796–799. [19] T.L. Christiansen, M.A.J. Somers, Defect Diffusion Forum 297–301 (2010) 1408–1413. [20] M. Bornert, T. Bretheau, P. Gilormini, Homogénéisation en Mécanique des Matériaux, vol. 2, Hermès, Paris, 2001. [21] I. Kaur, Y. Mishin, W. Gust, Fundamentals of Grain and Interphase Boundary Diffusion, Wiley, Chichester, 1995. [22] T. Roland, D. Retraint, K. Lu, J. Lu, Mater. Sci. Forum 524–525 (2006) 717–722. [23] E. Saatdjian, Phénomènes de transport et leurs résolutions numériques, Polytechnica, 1998. [24] V.H. Garcia, P.M. Mors, C. Scherer, Acta Mater. 48 (2000) 1201–1206. [25] H. Wei, X. Sun, Q. Zhen, G. Hou, H. Guan, Z. Hu, Surf. Coat. Technol. 182 (2004) 112–116. [26] J.M. Brossard, B. Panicaud, J. Balmain, G. Bonnet, Acta Mater. 55 (2007) 6586–6595. [27] T. Moskalioviene, A. Galdikas, J.P. Rivière, L. Pichon, Surf. Coat. Technol. 205 (2011) 3301–3306. [28] I.V. Belova, G.E. Murch, J. Phys. Chem. Solids 64 (2003) 873–878. [29] J.R. Kalnin, E.A. Kotomin, J. Maier, J. Phys. Chem. Solids 63 (2002) 449–456. [30] Y. Sun, T. Bell, Mat. Sci. Eng. A 224 (1997) 33–47.