Finite element simulation of diffusion into polycrystalline materials

Available online at www.sciencedirect.com

Solid State Sciences 10 (2008) 754e760 www.elsevier.com/locate/ssscie

Finite element simulation of diffusion into polycrystalline materials D. Gryaznov a,c, J. Fleig a,b, J. Maier a,* b

a Max Planck Institute for Solid State Research, Heisenbergstrasse1, D-70569 Stuttgart, Germany Vienna University of Technology, Institute of Chemical Technologies and Analytics, Getreidemarkt 9/164EC, A e 1060 Vienna, Austria c Institute for Solid State Physics, University of Latvia, Kengaraga 8, LV-1063 Riga, Latvia

Received 23 November 2007; received in revised form 27 March 2008; accepted 30 March 2008 Available online 18 April 2008

Abstract Diffusion in polycrystalline materials is investigated by means of numerical finite element simulations for constant source conditions. The grain boundaries are assumed to provide fast diffusion paths. Main emphasis is put on situations that typically occur for nanocrystals, viz. on situations in which (i) the diffusion length is significant compared with grain size, (ii) the influence of boundaries that are parallel to the surface become important in addition to the perpendicular ones. Furthermore, we treat the influence of blocking space charge layers sandwiching the core pathways and thus channeling grain boundary diffusion. Ó 2008 Elsevier Masson SAS. All rights reserved. Keywords: Diffusion; Grain boundaries; Nanocrystals; Finite elements; Space charges

1. Introduction Diffusion in polycrystalline materials is of crucial importance in solid state science and technology and frequently influences the preparation and characteristics of devices such as sensors and fuel cells. In many cases the grain boundary (GB) diffusion is decisive for the kinetics as the GB diffusion coefficient may significantly differ in comparison to that of the bulk and thus influences, for example, corrosion in metals and stoichiometry changes in oxides. In order to determine diffusion coefficients [1] from measured diffusion profiles the solutions of Fisher’s model [2] and in particular Le Claire’s relation [3] are commonly used to analyze experimental data. They assume polycrystalline microstructures with large grains and diffusion lengths in the grain (Lg) being much smaller than the average grain size (d ) [4]. The treatment of nanocrystalline materials, however, requires modifications of classical grain boundary diffusion models [5,6]. Several modifications of conventional models are discussed in the literature: Mishin [7], for

* Corresponding author. Tel.: þ49 711 689 1720; fax: þ49 711 689 1722. E-mail address: [email protected] (J. Maier). 1293-2558/$ - see front matter Ó 2008 Elsevier Masson SAS. All rights reserved. doi:10.1016/j.solidstatesciences.2008.03.030

example, analyzed the influence of the inclination of GBs to a free surface on the diffusion characteristics. He proposed a modified expression for the triple product sdDgb from which the GB diffusivity Dgb is to be extracted (s and d denote segregation factor and GB thickness). Recently, also some attempts were made to find analytical solutions to diffusion equations, taking into account realistic microstructures [8]. However, these analytical solutions are restricted to the socalled type-C diffusion kinetics, i.e. they are restricted to the initial stages in which only the grain boundaries become filled by the diffusion into the core of grain boundaries [9]. Many findings on GB diffusion are based on the simulations of spherical grain patterns, including the observation of a linear dependence of ln Cav and y6/5 by Levine and MacCallum [10], where Cav is the average concentration measured as a function of penetration y. An improved evaluation procedure of such profiles was recently proposed based on numerical calculations [5]. While the study of samples with essentially blocking grain boundaries is of great interest, too (see e.g. [11]), experimental and theoretical researches have traditionally focused on grain boundaries that provide fast transport paths which are the usual case in metallurgy. In ceramic systems space charges are often more resistive than the grain interior and thus play

D. Gryaznov et al. / Solid State Sciences 10 (2008) 754e760

an additional role [11]. Of particular interest in this context is the situation in which the fast core pathway is sandwiched between blocking space charge zones. We will examine such latter cases under conditions that are in particular more realistic in view of the interest in nanocrystalline matter. There, interaction of the grain boundaries in terms of their impact on the profile becomes soon important, and also the high number of boundaries parallel to the surface matters. 2. Grain boundary diffusion models and realization using the finite element method 2.1. Finite element representation of Fisher’s system Fisher’s system, i.e. 8  2  vCg ðx; y; tÞ v Cg ðx; y; tÞ v2 Cg ðx; y; tÞ d > > > ¼ D þ ; if jxj  g < vt vx2 vy2 2  vCgb ðy; tÞ v2 Cgb ðy; tÞ 2Dg vCg ðx; y; tÞ d > > > ¼ Dgb þ ; if jxj ¼ :  2 vt vy vx 2 d jxj¼d 2

ð1Þ was originally proposed for an isolated boundary model [4], in which a single GB in a semi-infinite sample is perpendicular to the surface. The x- and y-coordinates refer to the directions parallel and perpendicular to the surface. The isolated boundary model is symmetrical with respect to x ¼ 0 and reduces the problem to a two dimensional (2D) representation of a GB. In Eq. (1) d is the GB thickness, which was fixed to 0.5 nm for all calculations shown here, Cg(x,y,t) and Cgb( y,t) are the grain and GB concentration, respectively. Being based on Fick’s second law [12] for bulk and GB diffusion, Fisher’s system consists of a set of two differential equations with constant diffusion coefficients in grains (Dg) and along GBs (Dgb). Fick’s second law for the GB is included in a simplified form taking into account continuity conditions [4] and corresponding to diffusion along a line, i.e. the finite GB thickness and lateral concentration variations within the GB are neglected. In this model also laterally averaged concentrations (Cav) do not include contributions of the grain boundary but only average concentrations in the grain (see discussion below). In our simulations the value for Dg, viz. 2.95  104 (nm)2/s, is taken from an experimental study of diffusion in undoped nanocrystalline ZrO2 at 500  C [13]. Even though one can find excellent books on the application of the finite element method (FEM) to diffusion and heat conduction problems (see, for example, [14]), however, to the best of our knowledge only a single paper deals with the application of finite element method (FEM) to GB diffusion [15]. There Fisher’s model was also used and the diffusion profiles calculated for the model of square grains at long diffusion times t were compared to those for an isolated boundary. In the present study, we employ the finite element program FLUX-EXPERT (Simulog, France) to investigate the further aspects of grain boundary diffusion. Calculations of 2D distribution of Cg [16] require the definition of the expression

ZZ

755

Z   vCg ðx; y; tÞ d dxdy þ ai Dgb grad Cg ðx; y; tÞ grad½ai dy vt 2 dG Z d vCgb ðy; tÞ dy ai þ 2 vt dG ZZ   þ Dg grad Cg ðx; y; tÞ grad½ai dxdy ¼ 0 ð2Þ

ai G

G

which reflects the continuity condition at the grain boundarye grain interface. The symbol ai denotes the projection polynomial. This expression represents a weak formulation of Fick’s second law and includes the projection of the differential equation on the polynomials according to Galerkin’s projective method [17]. In Eq. (2) the second and the fourth integral represent the conduction parts, whereas the other two integrals are responsible for the time variation. The GB conduction part (the second integral) is obtained from Gauss’s theorem [18] and the GB diffusion equation given by Eq. (1), [19]. 2.2. Modification of Fisher’s system for space charge diffusion problems In ionic materials space charge layers (SCLs) are often found adjacent to the GB core with properties heavily differing from the grain interior as well as from the grain boundary core; in many experimental situations the SCL is depleted of mobile charge carriers (blocking effect) [11]. On the basis of the grain boundary model described by Jamniket al. [20], we thus extended Fisher’s system in the following way (Fig. 1): first, the GB (or the GB core in the SCL problem) is again considered to be very thin with negligible lateral concentration variations. Diffusion in the SCL is described by the diffusion coefficient Dscl and the concentration Cscl(x,y,t). The assumption of a constant (effective) diffusion coefficient in the entire space charge region is a severe simplification but still allows us to consider the basic effects of blocking space charges at grain boundaries. A flux of particles from the GB core into the SCL occurs, because Dgb is supposed to be larger than 0

Surface

x Dg Dscl cg

Dgb

Dg Dscl

δ δscl

δscl

cg

Cgb cscl

cscl

y Fig. 1. Schematic representation of the grain boundary geometry used in the case of calculations with SCLs.

D. Gryaznov et al. / Solid State Sciences 10 (2008) 754e760

756

Dscl and Dg. The combination with continuity conditions leads to a system of three equations:

2

ð3Þ

while Eq. (3) presents a system of well-defined differential equations adapted to our problem, and while also Eq. (2) can be used for FEM realization of our problem (just replacing Dg in the fourth integral by Dscl for the integration in the SCL), the number of mesh elements in the regions required to accurately integrate Eq. (3) can be very different depending on the corresponding diffusivities and diffusion time (t). The integration of Fisher’s system (Eq. (1)) or its modified form for space charge problems (Eq. (3)) provides a 2D distribution of concentrations within the grain and space charges. Since measured diffusion profiles typically represent average concentrations also the calculated concentrations are averaged according to 1 Cav ðy; tÞ ¼ xmax  xmin

Zxmax Cðx; y; tÞdx;

ð4Þ

xmin

where xmax  xmin is the width of the sample used with the SCLs being included in the integration. In accordance with Fisher’s model the concentration in the GB core is still excluded. This allows direct comparison with established models. Note, all concentration variations discussed below suppose constant source conditions and normalization with respect to the surface concentration. 3. Results and discussion 3.1. Accuracy of finite element calculations The finite element results obtained by using FLUX-EXPERT can easily be compared with an integration of Whipple’s solution [21] for an isolated GB perpendicular to the free surface (obtained by using MATLAB [22]). The latter is exact in the sense that it is an analytical solution of Fisher’s system. The two results are compared in Fig. 2. One can see that they coincide within a negligible error until the influence of the zeroflux boundary condition (bottom end of the sample) becomes significant (the corresponding derivative goes to zero in the FLUX-EXPERT result). This check demonstrates that the model applied in FLUX-EXPERT can indeed be used to numerically integrate Fisher’s system. Such tests were particularly important because diffusion was also studied under the extreme conditions of short and long t, covering both B- (initial diffusion from surface and grain boundaries into the grain interior) and A-regime (final stage of diffusion) [9]. In particular short times require very dense meshes and the number of (triangular) mesh elements can be as large as 500 000 [19].

MATLAB FLUX-EXPERT -0.003

∂lnCav /∂y6/5, nm-6/5

 2  8 vCg ðx; y; tÞ v Cg ðx; y; tÞ v2 Cg ðx; y; tÞ d > > ¼ D ; if jxj  þ dscl þ > g > > 2 vt vx2 vy2 > 2  > < vC ðx; y; tÞ v Cscl ðx; y; tÞ v2 Cscl ðx; y; tÞ d d scl ¼ Dscl ; if  jxj þdscl þ 2 2 vt vx vy 2 2 > >  > > > vCgb ðy; tÞ v2 Cgb ðy; tÞ 2Dscl vCscl ðx; y; tÞ d > > ¼ Dgb þ :  d ; if jxj ¼ 2 vt vy2 vx d jxj¼

0.000

-0.006

-0.009

-0.012

-0.015 0

300

600

900

1200

1500

1800

y6/5 , nm6/5 Fig. 2. Comparison of the diffusion profile calculated by the FEM and the profile according to Whipple’s solution integrated in the program MATLAB. The diffusion profiles were obtained for D ¼ 2.2  104 at t ¼ 2000 s (Diffusion length in the bulk z 0.8 nm).

According to Fig. 2 the FEM results may be influenced by the presence of the ‘‘reflecting boundary condition’’ at the bottom of the sample and the sample length (or diffusion time) has thus to be chosen properly. Shewmon [23] discussed that for an accurate experimental determination of the diffusivity (D) a profile with a concentration variation by at least one order of magnitude should be used, which corresponds to a minimum profile length of y z 3ODt. The same criterion was found to be applicable in the present study. In finite systems, such as thin films, the reflecting boundary is a necessary condition [24] but it can severely affect the analysis and a fitting of the profiles to a straight line according to standard analysis procedures can lead to an overestimation of Dgb. Particularly, in the case of blocking space charges very dense meshes are necessary and effects of the reflecting boundary become even more pronounced due to the great difference between Dgb and Dscl. In Fig. 3 a part of the mesh used to simulate diffusion including space charge regions is shown: The number of triangle elements in the SCL is about 87% of the whole number of the elements in the sample. 3.2. Comparison of the geometrical models of square grains and parallel boundaries The model of square grains used here to analyze the diffusion behavior in realistic microstructures represents a 2D pattern with GBs perpendicular or parallel to the diffusion direction (Fig. 4). Also a zero flux condition is used at the sample’s bottom as boundary condition. The typical grain sizes d used in the simulation were 10, 25, 50 and 100 nm, and, consequently, 5e50 serial grain layers are considered. As a second geometry the model of parallel boundaries (without perpendicular ones) is used and there d simply means the distance between two neighboring GBs. If not indicated otherwise, the typical diffusivity ratio,

D. Gryaznov et al. / Solid State Sciences 10 (2008) 754e760

SCL

Bulk

D¼

Fig. 3. Part of the mesh used for the simulations of diffusion including effects due to SCLs. The GB is located on the right hand side.

x

Free surface

757

Dgb Dg

is 2.2  104. All ratios used allowed for sufficient accuracy and convergence in the FEM simulations. It is shown in Fig. 5 that under these conditions (type B regime) the GB parts of the profiles for the parallel boundary model exhibit the same slopes independent of d, whereas the so-called bulk part (close to the surface and ideally given by a complementary error function) is affected by GBs. Only for larger d the bulk parts of the profiles are close to the complementary error-function solution (solid black curve in Fig. 5). The strongest deviation occurs for very small d suggesting that one may get serious errors in determining Dg from the corresponding profile part, especially at short times. When perpendicular GBs come into play, additional features can be observed. In Fig. 6 profiles calculated for the model of square grains are compared to those for parallel boundaries. The perpendicular GBs of square grains lead to spikes which reflect the particle flux from these perpendicular grain boundaries into the grains. However, not only spikes modify the profile but the average slope of the entire concentration profile (neglecting these spikes) also changes. Hence such a slope is no longer an accurate measure of Dgb. It is seen that the discrepancy between the two models is more pronounced for smaller d while the slope for square grains increases with decreasing d. The lower bulk concentration in the square grains (except spikes) can be understood when considering the perpendicular grain boundaries as a kind of leakage path for the diffusant in the GB network. Thus the concentration in the parallel GBs of the square grains is smaller compared to the case of parallel boundaries solely and accordingly also the bulk concentration in square grains is smaller (except close to the perpendicular boundaries, i.e. spikes).

y Fig. 4. Concentration distribution obtained by using the FEM for the model of square grains (d ¼ 10 nm). Bright colors reflect enhanced concentrations. In this case the conditions even satisfy the requirements of type-C diffusion kinetics. Only a part of the sample is shown.

Fig. 5. Variation of ln Cav with y6/5 calculated for the model of parallel boundaries for different distances between the boundaries (d ) and at t ¼ 2000 s. The black curve corresponds to pure bulk diffusion given by a complementary error function.

758

D. Gryaznov et al. / Solid State Sciences 10 (2008) 754e760

a

b

Fig. 7. Derivatives vln Cav/vy6/5 calculated for square grains with d ¼ 25 and 100 nm (only the grain next to the surface is considered, D ¼ 2.2  104 and t ¼ 2000 s).

the concentrations in both models are almost the same after long t. However, in that case the A-regime is already reached. In this regime other ways to determine the diffusion coefficient are discussed [4,25]. Interestingly the grain boundary volume fraction seems to be a problematic measure to weight Dg and Dgb in such analysis models since geometries with different fractions (square grains and parallel boundaries) can lead to

a x, nm

b 25

0

x, nm 475 500

Fig. 6. Comparison of the diffusion profiles calculated for the models of square grains and of parallel boundaries for two different d (D ¼ 2.2  104, t ¼ 2000 s).

525 550 575

In Ref. [5] we suggested to determine the maximum of the derivative of the profile in a ln Cavey6/5 plot in order to obtain an accurate grain boundary diffusion coefficient in the case of isolated grain boundaries (Whipple’s solution). Analyzing the derivatives of the profiles in the grain closest to the surface (Fig. 7), however, reveals that a maximum cannot be reached even for ‘‘large’’ grains i.e. d ¼ 100 nm. In the case of parallel boundaries (similar to an isolated boundary) the position of the 6/5 maximum is found for y6/5 at t ¼ 2000 s for max z 300 nm 4 D ¼ 2.2  10 while the positions of perpendicular GBs are at 15.85, 47.59, 251.19 nm6/5 for d ¼ 10, 25, 100 nm, respectively. Hence, methods for determining Dgb based on Le Claire’s procedure should further be improved to account for the GBs orientations. Moreover, the effect of perpendicular GBs is reflected in a strong decrease of the derivative close to the spikes. The difference between the profile shapes of the two models as well as the spikes vanish with increasing t, suggesting the profiles to coincide when t becomes sufficiently large. In the concentration distributions shown in Fig. 8 it is seen that

600 625 650

1 0.9375 0.875 0.8125 0.75 0.6875 0.625 0.5625 0.5 0.4375 0.375 0.3125 0.25 0.1875 0.125 0.0625 0.205E-21

25

0 475 500 525 550 575 600 625

1 0.9375 0.875 0.8125 0.75 0.6875 0.625 0.5625 0.5 0.4375 0.375 0.3125 0.25 0.1875 0.125 0.0625 0.925E-28

650

675

675

700

700

725

725

750

750

775

775

800

800

y, nm

y, nm

Fig. 8. Section of the concentration distribution in some distance from the surface for t ¼ 3  106 s (type A-kinetics): (a) parallel boundaries, (b) square grains. White lines represent the perpendicular GBs, parallel boundaries in both distributions are at x ¼ 0.0 and 25.0 nm.

D. Gryaznov et al. / Solid State Sciences 10 (2008) 754e760

the same effective diffusion coefficient. More information on this topic is given in [19] and will be also discussed in a forthcoming paper.

layers in between which then the fast core path is embedded. In Ref. [6] we have shown on the basis of the FEM calculations for an isolated boundary model that for an increasing ratio,

3.3. Space charge effects

a

1.5

1.0

0.5

0.0 0.0

x, nm

1 0.9375 0.875

0.25

0.75 0.6875 0.625 0.5625

SCL

Grain

0.8125

0.50

0.5 0.4375 0.375

0.75

0.3125 0.25 0.1875 0.125 0.0625 0.154E-43

1.00

y,nm

b

6.0

4.0

2.0

0.0 0.0

x, nm SCL 2.0

4.0

6.0

1 0.9375 0.875 0.8125 0.75 0.6875 0.625 0.5625 0.5 0.4375 0.375 0.3125 0.25 0.1875 0.125 0.0625 0.154E-43 0

L¼

Dg Dscl

the slope of the profile ln Cav ¼ f( y6/5) decreases while for decreasing the opposite trend is observed. Accordingly, one has to consider the two cases of L < 1 and L > 1. The difference between the two situations is demonstrated in Fig. 9. If L > 1 (blocking boundaries), a drastically reduced concentration within the SCL results, what strongly affects the laterally averaged diffusion profile obtained under such conditions. One also can see that the isoconcentration lines within the SCL in Fig. 9a are almost parallel to the GB core, i.e. the angle between them is very small e a situation typical for short diffusion lengths in this region. According to this, the SCL is not only filled from the surface and the grain boundary core but also by a particle flux from the grain interior into the SCL. If on the other hand, L < 1 (Fig. 9b), then there is a flux from the SCL into the grain. The perpendicular GB present in the model of square grains further affects the slope and shape of the diffusion profile, especially at short t and/or large ratios between the diffusivities. In Fig. 10 the profile calculated from the model of square grains with SCL is compared with those obtained from the models of parallel boundaries and square grains without SCLs; in all cases L ¼ 102 and D ¼ 102 are used. Owing to the smaller D-value, spikes are less pronounced than in Fig. 6. Accordingly, the profile for the square grains without SCLs is close to the profile for parallel boundaries without SCLs. However, the slope of the profile for square grains with SCL is much smaller. This can be understood when taking into account that for identical surface concentrations the blocking character of the SCLs leads to a reduced ‘‘leakage’’ of diffusants from the GB cores into the grains. This results in 0 square grains with SCL -5

parallel boundaries without SCL

-10

lnCav

For studying the influence of space charges we approximate the space charge zones by homogeneous thin layers with modified but constant diffusion coefficient Dscl (For the connection between an averaged value with the concentration profile in the case of transport parallel and perpendicular to the interface see Ref. [26]). We will concentrate on the case of depletion

759

square grains without SCL -15

8.0 -20

y, nm

-25 0

Fig. 9. Section of the concentration distribution for the model of a single (isolated) GB exactly at x ¼ 0 nm with adjacent SCLs of thickness 1 nm. The free surface is located at y ¼ 0 nm. Parameters: (a) L ¼ 103, D ¼ 102, t ¼ 4700 s and (b) L ¼ 101, D ¼ 102, t ¼ 4700 s. Please note the different length scales in (a) and (b).

100

200

300

400

500

y6/5, nm6/5 Fig. 10. Variation of ln Cav with y6/5 calculated for a sample with square grains and GBs surrounded by blocking SCLs and, for comparison, for square grains and parallel boundaries without SCLs for L ¼ 102 and D ¼ 102 at t ¼ 8200 s.

D. Gryaznov et al. / Solid State Sciences 10 (2008) 754e760

760 0

25 x, nm

1 0.937507

25

0.875013 0.81252 0.750027 0.687534 0.62504 0.562547

50

0.500054 0.437561 0.375067 0.312574 0.250081

75

0.187588 0.125094 0.062601

proved particularly helpful when addressing nanocrystalline materials. This sensible application requires high mesh densities and a careful analysis of the influence of boundary conditions. On the basis of simulations by FEM it was observed that GBs perpendicular to the diffusion flux play an important role when the GB diffusion length is larger than the average grain size. Only for very long times, the diffusion profiles largely coincide for the models of parallel boundaries and square grains. For blocking space charge layers being present in ionic materials the conventional determination of Dgb from the diffusion profiles can be rather problematic since the channeling of core transport severely modifies the concentration distribution. Moreover, blocking SCLs lead to abrupt steps of the concentration at perpendicular GBs, and extremely long times are required to fill all the grains of ionic materials with blocking space charge layers at grain boundaries. At any rate a neglect of such phenomena would lead to serious misinterpretations in terms of kinetic constants.

0.000108

References

100 y, nm

Fig. 11. Section of the concentration distribution in the case of a sample with square grains and SCLs at the GBs for t ¼ 3  106 s and L ¼ 103. The GBs are located at x ¼ 0.0 and 25.0 nm, y ¼ 25, 50, 75 and 100 nm.

a kind of channeling of the diffusants and their penetration depth in the GB core becomes larger. Accordingly, in some depth the (low) flux from the GB core into the grains becomes higher than in samples without SCL which is simply an effect of the increased boundary concentration in the core. The concentration of diffusants close to the surface, however, indicates that despite the deep penetration, the total amount of diffusant entering the sample is decreased by the blocking space charge. For long diffusion times another feature of blocking SCLs becomes visible (Fig. 11): While the top grain layer becomes completely filled by bulk diffusion from the surface other grains are characterized by much smaller concentrations. Accordingly, a sharp step between first and second grain layer appears in the profile. This is due to the fact that the perpendicular GBs with SCLs do hardly admit any particle flux from the first to the second layer; in addition the leakage of particles from the fast GB cores into the grains is minimized due to blocking SCLs. 4. Conclusions The finite element method was used to simulate diffusion in polycrystals under various conditions of relevance. It

[1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26]

See e.g. H. Schmalzried, Solid State Reactions, VCH, Weinheim, 1981 J.C. Fisher, J. Appl. Phys. 22 (1) (1951) 74. A.D. Le Claire, Brit. J. Appl. Phys. 14 (1963) 351. Kaur, Y. Mishin, W. Gust, Fundamentals of Grain and Interface Boundary Diffusion, Wiley, Chichester West Sussex, 1995. D Gryaznov, J Fleig, J Maier, J. Appl. Phys. 103 (2008) 063717. D. Gryaznov, J. Fleig, J. Maier, Solid State Ionics 177 (2006) 1583. Yu. M. Mishin, Phys. Status Solidi A 133 (1992) 259. K. Bedu-Amissah, J.M. Rickman, H.M. Chan, M.P. Harmer, J. Appl. Phys 98 (2005) 063511. L.G. Harrison, Trans. Faraday Soc. 57 (1961) 1191. H.S. Levine, C.J. MacCallum, J. Appl. Phys. 31 (3) (1960) 595. J. Maier, Physical Chemistry of Ionic Materials. Ions and Electrons in Solids, John Wiley & Sons, Chichester, 2004. A. Fick, Philos. Mag. J. Sci. 4 (10) (1855) 30. U. Brossman, R. Wu¨rschum, U. So¨dervall, H.-E. Schaefer, J. Appl. Phys. 85 (11) (1999) 7646. G. Comini, S.D. Guidice, C. Nonino, Finite Element Method in Heat Conduction, Taylor & Francis, Washington, 1994. Z. Knesl, J. Jinoch, Z. Bilek, Czech. J. Phys. B 24 (1974) 347. FLUX-EXPERTÒ (Simulog), User Manual, Version 3.1, 2000. Handbook on Finite Element Method and FLUX-EXPERT (DT2i), Introduction to the Finite Element Method, Version 1.6, 1992. G.A. Korn, T.M. Korn, Mathematical Handbook for Scientists and Engineers, McGraw-Hill, New York, 1968. D. Grjaznovs, PhD thesis, University of Stuttgart, 2006. J. Jamnik, J. Maier, S. Pejovnik, Solid State Ionics 80 (1995) 19. R.T.P. Whipple, Philos. Mag. 45 (1954) 1225. MatLabÒ, The Language of Technical Computing, Mathematics (The Mathworks, Inc), Version 7, 2004. P.G. Shewmon, Diffusion in Solids, McGraw-Hill, New York, 1963. G.H. Gilmer, H.H. Farrel, J. Appl. Phys. 47 (10) (1976) 4373. I.V. Belova, G.E. Murch, Philos. Mag. 81 (10) (2001) 2447. J. Maier, Ber. Bunsenges. Phys. Chem. 90 (1986) 26.