- Email: [email protected]
Numerical modelling of hydrogen desorption from cylindrical surface
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 6 ( 2 0 1 1 ) 1 2 3 9 e1 2 4 7
Available at www.sciencedirect.com
journal homepage: www.elsevier.com/locate/he
Numerical modelling of hydrogen desorption from cylindrical surface Natalia I. Rodchenkova*, Yury V. Zaika Institute of Applied Mathematical Research, Karelian Research Centre of the Russian Academy of Sciences, Pushkinskaya 11, 185910 Petrozavodsk, Russian Federation
article info
abstract
Article history:
We consider degassing of cylindrical sample containing dissolved hydrogen. The experi-
Received 25 June 2010
ment is made by the thermodesorption method. In the corresponding boundary-value
Accepted 25 June 2010
problem with nonlinear dynamic boundary conditions physicalechemical processes in the
Available online 17 September 2010
bulk and on the metal surface are taken into account: diffusion, desorption, solution, and capture by defects. Computational algorithm for desorption flux modelling is developed on
Keywords: Hydrogen thermodesorption
the basis of difference approximations. The results of numerical modelling are presented. ª 2010 Professor T. Nejat Veziroglu. Published by Elsevier Ltd. All rights reserved.
Mathematical modelling
1.
Introduction
A high concentration of hydrogen in metal leads to hydrogen embrittlement [1,2]. Common metallurgical concentrations of dissolved hydrogen range from 0.1 to 100 ppm. The authors of paper [3] have developed the hydrogen analyzer AV-1. The device allows to measure the concentration of hydrogen atoms (H) in a solid sample using industrial laboratory facilities. A cylindrical sample is placed into the quartz vacuum extractor. The extractor is exposed to a given temperature in the furnace. The contact between the sample and the extractor walls is point-like, the thermal conductivity of quartz is negligible, wherefore the heat transfer is determined by radiation. As the sample is heated, atomic hydrogen diffuses inside the bulk and is desorbed from the surface in the molecular form. The extraction curve (measured by a mass-spectrometric hydrogen analyzer) is recorded. The curve is used for further processing (in particular, the kinetic parameters of models are estimated). The dependence of the desorption flux on the temperature at uniform heating (the TDS-spectrum) usually has several peaks. According to [4], in addition to the diffusion, other surface processes may affect
the desorption of hydrogen (e. g. the capture of hydrogen atoms by cracks, microcavities; inclusion of hydrogen into hydride phases). The experimental approaches to the problems of hydrogen materials science, especially in the case of heavy isotopes deuterium and tritium, are rather complicated and expensive. Therefore, the role of mathematical modelling in hydrogen materials science is very important. The paper is devoted to the development of mathematical tools for reliable experimental measurements of the hydrogen content in various materials. In this work, we will use parameters suitable for aluminium and some of its alloys when numerically modelling the desorption flux.
2.
Mathematical model
2.1.
Heating equation
The sample has the form of a cylinder with typical dimensions (GOST 21132.1-98): radius of the base is L ¼ 4 103 m, height is H ¼ 2102 m. Concentration of diffusing atomic hydrogen at the initial time instant cð0; r; zÞ ¼ c is constant (it forms during
* Corresponding author. Tel.: þ7 8142766312; fax: þ7 8142766313. E-mail addresses: [email protected] (N.I. Rodchenkova), [email protected] (Yu.V. Zaika). 0360-3199/$ e see front matter ª 2010 Professor T. Nejat Veziroglu. Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.ijhydene.2010.06.121
1240
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 6 ( 2 0 1 1 ) 1 2 3 9 e1 2 4 7
Nomenclature L, H S, V t ts T Te tcrit, Tcrit r s ~c l g
radius, height of the cylinder (m) surface area and capacity (m2, m3) time (s) time instant when surface H is exhausted (s) temperature (K) temperature of extractor wall (K) critical time and temperature (s, K) bulk density of the metal (kg/m3) Stefan-Boltzmann constant (J/s/m2/K4) specific heat capacity (J/kg/K) thermal conductivity (J/s/m/K) quick dissolution coefficient (1/m)
the material production). When necessary, it is possible to take account of the decrease in H concentration in the near-surface layer (for instance, as the result of mechanical and thermal pretreatment) without fundamental changes in the numerical algorithm. If the sample heating is uniform (sufficiently slow, T ¼ T(t), ½T ¼ K), the temperature dynamics can be described by the differential equation of the heat balance [5]: dT sS 7 105 ðT þ 64:3Þ T4e T4 ; ¼ dt ~crV
T0 ¼ Tð0Þ ¼ 293 K:
Here, Te ¼ const is the temperature of the extractor wall, S ¼ 2pLH, V ¼ pL2H are the surface area and the cylinder capacity, r ¼ 2.71 103 kg m3 is the bulk density, s ¼ 5.67 108 J s1 m2 K4 is the Stefan-Boltzmann constant, ~c ¼ 1:15 103 J kg1 K1 is the specific heat capacity. Since the assumption about uniform heating does not have adequate accuracy for all values of Te, we consider an alternative distributed model. We suppose that in the case of tubular geometry of the extractor heating mainly proceeds through the lateral surface. Thus, we speak of a lower estimate of the heating of the centre of the sample. The considered equation from [5] is majorized from above. The radiosymmetrical model is: T(0, r) ¼ T0 < Te, r˛½0; L,
J hydrogen desorption flux (1/s) k thermal diffusivity (m2/s) heat absorption coefficient 3t a emissivity (W/m2/K4) R universal gas constant (J/mol/K) r, z cylindrical coordinates (m) c H concentration in metal (1/m3) w H concentration in defects (1/m3) g H concentration in hydride phase (1/m3) c; w; g initial concentrations (1/m3) surface concentrations (1/m2) qi defect coefficients (1/s) ai E activation energy (J/mol) D, b diffusion and desorption (m2/s)
2.2.
Diffusion model with defects
Consider the boundary-value problem of thermodesorption for the cylinder taking into consideration diffusion in the bulk, H capture by defects, penetration from the bulk to the surface and desorption: 2 vc v c 1 vc v2 c þ $ þ 2 a1 ½1 Wcðt; r; zÞ ¼DðTÞ 2 vt vr r vr vz (1) þ a2 wðt; r; zÞ þ a3 gðtÞ; vw w ; ¼ a1 ðTÞ½1 Wc a2 ðTÞw; Wh vt wmax (2) r˛ð0; LÞ; z˛ð0; HÞ; t˛ð0; t Þ;
2 vT v T 1 vT vT vT þ ¼ 0; l ¼ a T4e T4 ðt; LÞ : ¼k $ ; vt vr2 r vr vr r¼þ0 vr L
dg ¼a3 gðtÞ; gðtcrit Þ¼g0gðtÞ¼gexpfðtcrit tÞa3 g; a3 ¼0; T
Here, l is the heat conduction coefficient (for aluminium in the range of T˛½300; 800 lz236 J s1 m1 K1 ); k ¼ l ð~crÞ1 is the thermal diffusivity; the finishing time of computing t* is determined by the stationary distribution Tðt; 0ÞzTe , t > t*; a ¼ s3t, 3t ¼ 7 105(T þ 64.3) is the heat absorption coefficient. The auxiliary problem of numerical modelling of heating is the following: it is necessary to quickly estimate how much the distribution T(t, r) differs from the uniform heating T(t) under given Te, L, H, and thermal and physical characteristics of the material. For instance, when we use T0 ¼ 293 K, Te ¼ 773 K and the L, H stated above the assumption about uniform heating of aluminium sample is true: the difference T(t) T(t, 0) is within a hundredth of a degree. We obtain this result in the distributed model where heat absorption by butt ends is disregarded. Peak temperature is reached in 2.2 h, heating being practically uniform within an hour. Hereinafter, suppose that the sample heating is uniform.
Here, c(t, r, z) is the concentration of corpuscular (diffusion movable) hydrogen in the metal; w(t, r, z) is H concentration in the defects of the material physicalechemical structure, the capture is reversible (for instance, microcavities); wmax is the peak concentration of reversible capture; g(t) is H concentration in traps, which begin to release the hydrogen only when some critical temperature Tcrit is reached (it is typical for the inclusion of hydride phases); ai are the coefficients of absorption and desorption of H by traps (a3 > 0 only at T Tcrit); q1(t, z), q2(t, r), q3(t, r) are the surface concentrations (on the lateral surface of the cylinder and on butt ends); c is the initial concentration c(0,r,z) in the solid sample; g is the fit coefficient for concentrations of hydrogen atoms in the bulk and on the surface (quick dissolution coefficient); D, b are the diffusion and the desorption coefficients. Suppose that ai > 0 are constant within the considered temperature range ðT˛½300; 800Þ. Changes in the
1241
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 6 ( 2 0 1 1 ) 1 2 3 9 e1 2 4 7
time-dependent case of ai ðtÞhai ðTðtÞÞ are insignificant. For practical purposes we consider capture in a simple integrated form. Specifications of defects geometry and their distributions in the sample essentially complicate the model. To simplify the computations we consider one generalized trap with reversible capture and one generalized trap with decomposition. It is possible to consider several ðjÞ ðkÞ traps with their individual ai and Tcrit . We assume that diffusion and desorption coefficients depend on the temperature in an Arrhenius way: D(T ) ¼ D0exp{ED/[RT]}, b(T ) ¼ b0exp{Eb/[RT]}, ED,b are activation energies. We use a shorthand notation D(t) ¼ D(T(t)), b(t) ¼ b(T(t)). The time segment [0, t*] is determined by total degassing of the sample: JðtÞz0, t t*. The conditions cr(t, þ0, z) ¼ 0, cz(t, r, H/ 2) ¼ 0 are consequences of symmetry. A more precise model of dissolution on the surface (considering the lateral surface for clarity) has a form of fluxes balance: 1 kþ ðTÞcðt; L; zÞ 1 q1 ðt; zÞq1 max k ðTÞq1 ðt; zÞ 1 cðt; L; zÞcmax ¼ DðTÞcr ðt; L; zÞ: However, when diffusion is considerably slower than dissolution (the temperature is not very low) and concentrations are low, we obtain the condition of quick dissolution czgq, where g ¼ k/kþ. If the surface is isotropic (in the sense Ek zEkþ ), then the parameter g weakly depends on T. Below we use the notation Eg conditionally: it is not the activation energy (the difference Ek Ekþ may be negative). For traps which begin to release hydrogen only when certain critical temperature is reached (traps of the hydride phase type) we did the following: we took into account only the trap capacity and decomposition rate. Modelling of the dehydriding process is an independent intricate problem, which results in nonlinear boundary-value problems with free phase boundaries and with conditions of the Stefan type. The variant with bulk desorption is considered in Ref. [6], the variant with surface desorption is presented in Ref. [7]. Since the initial data are symmetrical we have q3 ¼ q2, and so below we construct the difference approximation only for a half of the cylinder ðz˛½H=2; HÞ with corresponding boundary conditions ðcz jH=2 ¼ 0; q_ 2 ¼ .Þ. For the defect with reversible capture (microcavity) the constant w is determined according to (2) ðv=vt ¼ 0Þ : a1 ðT0 Þ½1 w=wmax c a2 ðT0 Þw ¼ 0. For the trap of the hydride phase inclusion type the values of g ¼ const, Tcrit, a3 are set using information about specific chemistry of the hydride. The presence of derivatives q_ i (accumulation on the surface) corresponds to the possibility of H migration along the surface until H atoms form H2 molecules which are desorbed from the surface. Thus, the model aimed to analyze the problem of the dynamics of low common concentrations of hydrogen (without preliminary saturation in laboratory). The problem has an applied context, so the authors apply minimal mathbased environment to describe basic processes. Further detailed elaboration leads to an increase in the number of parameters. The inverse problem of the parameters estimation becomes complicated.
2.3.
Specification of the problem statement
The aim of the paper is to develop the difference scheme and computational algorithm for modelling of hydrogen desorption flux from cylindrical sample: Z JðtÞ ¼ 4pbðtÞ L
Z
H H=2
q21 ðt; zÞdz
þ
L 0
! rq22 ðt; rÞdr
:
The H2 molecule is desorbed from the surface, but we compute the desorption flux in atoms ([J] ¼ 1/s). The criterion of computation correctness is the material balance: Z pL2 Hðc þ w þ gÞ þ 4p L Z ¼ 4p
Z
H
Z þ 4p L
Z q1 ð0; zÞdz þ
H=2 L
dz H=2
H
! rq2 ð0; rÞdr
0
r½cðt; r; zÞ þ wðt; r; zÞ þ gðtÞdr
0 H
L
Z q1 ðt; zÞdz þ
H=2
L
! rq2 ðt; rÞdr
0
Z þ
t
JðsÞds: 0
Under uniform heating, it is convenient to consider the TDS-spectrum (the curve J ¼ J(T )) together with the flux dependence on t. TDS-spectrum usually has several peaks. It is believed the first peak corresponds to surface hydrogen. But one should have caution: while the surface hydrogen is desorbed, hydrogen atoms penetrate from the bulk to the surface. The problem of estimating the corresponding correction is topical. Let ts be the time instant when the surface hydrogen is exhausted. At the numerical modelling stage the time instant ts is determined from Z
ts
Z JðsÞdsz4p L
0
H H=2
Z q1 ð0; zÞdz þ
L
! rq2 ð0; rÞdr :
0
Numerical modelling allows to choose a region on the extraction curve which corresponds to the initial quantity of the surface hydrogen, to evaluate the value of the activation energies of diffusion and desorption, to estimate the parameters of reversible capture and the parameters of hydride phases decomposition.
3. Difference approximation of boundaryvalue problem 3.1.
Difference scheme construction
Follow the standard method [8], consider a spatial grid Uh ¼ ri ; zj : ri ¼ ihr ; i ¼ 0; 1; .; N1 ¼ ½L=hr ; zj ¼ jhz ;
j ¼ 0; 1; .; N2 ¼ ½ðH=2Þ=hz and a time grid us¼{tk ¼ ks, k ¼ 0, 1,., K ¼ [t*/s]}. Denote approximate values of the bulk concentration c(tk, ri, zj) by cki;j . Similarly, wki;j zwðtk ; ri ; zj Þ; gk zgðtk Þ; aks ¼ as ðtk Þ, Dk ¼ D(tk), where ðri ; zj Þ˛Uh , tk ˛us . For equation (1) examine the implicit difference scheme of the alternating directions method, called wedge-reed (PeacemaneRachford) scheme, while for equation (2) consider the scheme with weights. Transition from the kth layer to the (k þ 1)th layer is realized in two steps. At the first kþ1=2 are determined from the step intermediate values of ci;j equations set
1242
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 6 ( 2 0 1 1 ) 1 2 3 9 e1 2 4 7
kþ1=2 ci;j cki;j
0:5s
kþ1=2 kþ1=2 kþ1=2 ciþ1;j 2ci;j þ ci1;j
kþ1=2 kþ1=2 1 ciþ1;j ci1;j ¼ Dkþ1=2 þ $ h2r 2hr ri k k k ci;jþ1 2ci;j þ ci;j1 kþ1=2 kþ1=2 ~ þ Dk ci;j a1kþ1=2 1 W i;j h2z
!
kþ1=2 þ a2kþ1=2 wi;j þ a3 gkþ1=2 ;
kþ1=2 wi;j wki;j
0:5s
ð5Þ
h i k ¼ð1 sÞ ak1 1 Wi;j cki;j ak2 wki;j ~ kþ1=2 ckþ1=2 akþ1=2 wkþ1=2 : þ s akþ1=2 1 W 1
kþ1=2 kþ1=2 kþ1=2 kþ1=2 Ai ci1;j Gi;j ci;j þ Bi ciþ1;j þ Fki;j ¼ 0; k 0:
i;j
i;j
2
i;j
kþ1=2 ci;j
0:5s
kþ1=2 ciþ1;j
þ
kþ1=2 ci1;j
kþ1=2 ciþ1;j
kþ1=2 ci1;j
1 þ $ 2hr ri kþ1 kþ1 ckþ1 i;jþ1 2ci;j þ ci;j1 kþ1 ~ kþ1 ckþ1 1 W þ Dkþ1 a i;j 1 i;j h2z
¼Dkþ1=2
kþ1=2 2ci;j h2r
kþ1 kþ1 þ akþ1 ; 2 wi;j þ a3 g
kþ1=2
wkþ1 i;j wi;j 0:5s
ð6Þ
!
i;j
At r / þ0 we have cr =r ¼ ðcr ðt; r; zÞ cr ðt; 0; zÞÞ=rzcrr . Find the first coefficients from approximation of the equation
1 kþ1=2 ~ kþ1=2 ; bkþ1 ¼ Fk =4: 1þV cr jþ0 ¼ 0 : a1;j ¼ 1 k 4Dkþ1=2 1;j 1;j 1;j which is The immediate task is to find a value of cNkþ1=2 1 ;j necessary to launch the sweep method. Write down the approximation of the boundary condition (3) (r ¼ L) : (k 0),
. h
2 ð0:5sÞ ¼ 0:5 bk qk1j
2 i kþ1=2 Dk cr tk ; L; zj bkþ1=2 q1j Dkþ1=2 cr tkþ1=2 ; L; zj :
kþ1=2
q1j
qk1j
(9)
Sweep method with respect to the radius r
ð12Þ
In the boundary node with a precision of up to Oðh2r Þ we have m m 2hr cr tm ; L; zj zcm N1 2;j 4cN1 1;j þ 3cN1 ;j ; m ¼ k; k þ 1=2:
In the standard notations y_ ¼ f ðt; yÞ, ys ¼ ys1/2 þ {fs1/2 þ fs}s/4, the equation (9) is the symmetrical scheme for equation (2) normalized to wmax with the time-fixed function c ¼ c(ts1/2, r, z). We explain below the iterative procedure of more precise ~ s . For distinctness define s ¼ 1/2. The approxdefinition of W i;j imation order is Oðs2 þ h2r þ h2z Þ [8].
(13)
Concentration values on the kth layer are already known. For and cNkþ1=2 from the (k þ 1/2)th layer, substituting values cNkþ1=2 1 2;j 1 1;j expression (11), we have
1 h 3 þ aNkþ1=2 aNkþ1=2 cr tkþ1=2 ; L; zj z cNkþ1=2 þ bNkþ1=2 ;j 1;j 4 1 1 1 ;j 1 1;j 2hr i
: þ aNkþ1=2 4 bNkþ1=2 1 1;j 1 ;j Write down the approximation in a more compact way:
kþ1=2 kþ1=2 kþ1=2 kþ1=2 kþ1=2 Ah3 þ aN1 ;j aN1 1;j 4 ; BhbN1 1;j þ aN1 1;j 4 bN1 ;j ; þB: 2hr cr tkþ1=2 ; L; zj zAcNkþ1=2 1 ;j
ð14Þ
th
Consider a transition from the k layer to the (k þ 1/2) layer. kþ1=2 from the equation (6) and substitute it into (5). Express wi;j In notations Ai ¼ 1hr ð2ri Þ1 ; Bi ¼ 1þhr ð2ri Þ1 ; k ¼ 2h2r s1 ; h i k ak1 1Wi;j kþ1=2 kþ1=2 k ~ ; Vi;j ¼ 2þkD1 ¼ kþ1=2 ; Gi;j kþ1=2 1þ Vi;j þ4s1 a2 kþ1=2 1 ~ kþ1=2 1 W 2a1 i;j 4s ak2 wki;j kþ1=2 k ~ V ¼ ; N ¼ ; i;j i;j 4s1 þa2kþ1=2 a2kþ1=2 þ4s1
Dk 2 Fki;j ¼ ðhr =hz Þ cki;jþ1 2cki;j þcki;j1 Dkþ1=2 h i h i kþ1=2 k 1 k cki;j þh2r D1 þkDkþ1=2 1þ0:5sa2kþ1=2 Vi;j Ni;j þa3 gkþ1=2 ; kþ1=2 a2 at every fixed j ¼ 1, 2,., N2 1 obtain
. . kþ1=2 kþ1=2 kþ1=2 kþ1=2 ¼Bi1 Gi1;j Ai1 ai1;j ; bi;j ¼ Ai1 bi1;j þ Fki1;j
kþ1=2 kþ1=2 ; i ¼ 2; 3; .; N1 : Gi1;j Ai1 ai1;j
on the (k þ 1/2)th time layer for i ¼ 1 and using the condition
h i
kþ1=2 kþ1=2 kþ1=2 ci;j ¼ð1 sÞ a1kþ1=2 1 Wi;j a2kþ1=2 wi;j ~ kþ1 ckþ1 akþ1 wkþ1 : 1 W ð8Þ þ s akþ1 i;j 1 2 i;j i;j
th
kþ1=2
ai;j
ct ¼ Dð2crr þ czz Þ a1 ½1 Wc þ a2 w þ a3 g
Here, in order to have a possibility to use the sweep method on s the sth time layer (s ¼ k þ 1/2; k þ 1), the unknown Wi;j is replaced by its approximation from the linear equation with ~s respect to W
3.2.
(11)
Sweep coefficients are:
ð7Þ
h i n ~ s ¼ Ws1=2 þ s as1=2 1 Ws1=2 cs1=2 w1 as1=2 Ws1=2 W i;j 2 max i;j i;j i;j i;j 4 1 h i o s s ~ : ~ cs1=2 w1 as W þ as1 1 W i;j i;j i;j max 2
Values at the initial time instant (on zero layer) are known: c0i;j ¼ c ¼ const. According to the sweep method, search approximate values of the concentration in grid points on the (kþ1/2)th time layer in the form kþ1=2 kþ1=2 kþ1=2 kþ1=2 ci;j ¼ aiþ1;j ciþ1;j þ biþ1;j ; i ¼ 0; 1; .; N1 1; k 0:
kþ1=2 At the second step, using the resultant ci;j , find ckþ1 from i;j the system
ckþ1 i;j
(10)
Expressions (13), (14) substitute into (12), denote sNkþ1=2 ¼ y: 1 ;j 1 m qm 1j ¼ gm cN1 ;j ; m ¼ k; k þ 1=2; # " 2 bkþ1=2 2 Dkþ1=2 A 4 bk y þ G ¼ 0; Gh 2 ckN1 ;j y þ þ 2hr sgkþ1=2 g2kþ1=2 gk
3Dk 4 1 Dkþ1=2 B þ Dk ckN1 2;j 4ckN1 1;j : ckN1 ;j þ þ 2hr sgk 2hr
This is a square equation with respect to variable y. The >0 roots have different signs (s 1). We choose y ¼ sNkþ1=2 1 ;j taking into account the physical meaning of the process. Approximation order of the boundary condition is Oðh2r þ s2 Þ. It agrees with the approximation order of the scheme in the bulk.
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 6 ( 2 0 1 1 ) 1 2 3 9 e1 2 4 7
kþ1=2
Now find all values of ci;j at j ¼ 0 and j ¼ N2, i ¼ 0, 1,., N1. Using the symmetry condition cz jH=2 ¼ 0, we have kþ1=2 kþ1=2 kþ1=2 kþ1=2 ¼ ð4ci;1 ci;2 Þ=3. Values ci;N2 > 0 (s, h are small) are ci;0 unambiguously defined from the square equation which approximates the boundary condition (4) at z ¼ H: # " bkþ1=2 kþ1=2 2 3Dkþ1=2 kþ1=2 2 ci;N2 þ þ c 2hz sgkþ1=2 g2kþ1=2 i;N2 2 Dkþ1=2 kþ1=2 kþ1=2 þ ci;N2 2 4ci;N cki;N2 ¼ 0: 2 1 2hz sgk kþ1=2 , compute the concenWhen we know all values of ci;j tration of w:
h i kþ1=2 k wi;j ¼ 4s1 wki;j þ ak1 1 Wi;j cki;j ak2 wki;j 1 ~ kþ1=2 ckþ1=2 akþ1=2 þ 4s1 þ a1kþ1=2 1 W : 2 i;j i;j
Sweep method with respect to the variable z
3.3.
Since there is a singularity in cylindrical coordinates at r / þ0, we make a transition from the (k þ 1/2)th layer to the (k þ 1)th layer in two steps. First step is: i ¼ 1, r / þ0 (approximate cr =rzcrr ), we realize the sweep method for the equation ct ¼ D(2crr þ czz) a1[1 W]c þ a2w þ a3g. Second step is: i ¼ 2,., N1 1, r > 0, realize the sweep method for the equation ct ¼ (crr þ cr/r þ czz) a1[1 W]c þ a2w þ a3g. The equation for w remains the same. Equation (7) for i ¼ 1 now is kþ1=2 ckþ1 1;j c1;j
0:5 s
¼2Dkþ1=2 þ þ
kþ1=2 kþ1=2 kþ1=2 c2;j 2c1;j þ c0;j
h2r
kþ1 2ckþ1 1;j þ c1;j1 Dkþ1 2 hz kþ1 akþ1 w þ a gkþ1 : 3 2 1;j
ckþ1 1;jþ1
Difference approximation of equation (7) for i ¼ 2,., N1 1 is: kþ1=2 ckþ1 i;j ci;j
kþ1=2 kþ1=2 1 ciþ1;j ci1;j þ $ 2 0:5s hr 2hr ri kþ1 kþ1 ckþ1 2c þ c i;jþ1 i;j i;j1 ~ kþ1 ckþ1 1W þ Dkþ1 akþ1 i;j 1 i;j h2z
¼ Dkþ1=2
kþ1=2 kþ1=2 kþ1=2 ciþ1;j 2ci;j þ ci1;j
kþ1 kþ1 þ akþ1 : 2 wi;j þ a3 g
!
ð18Þ
Find wkþ1 from (8) (s ¼ 1/2) and substitute it into (18). In i;j notations h i kþ1=2 a1kþ1=2 1Wi;j kþ1 kþ1=2 2 1 kþ1 1 ~ ; Vi;j ¼ ; k ¼ 2hz s ; Gi;j ¼ 2þkDkþ1 1þ V i;j akþ1 þ4s1 2
~ kþ1 kþ1=2 1 W 2akþ1 1 i;j 4s1 a2kþ1=2 akþ1 2 wi;j kþ1 kþ1=2 ~ ¼ V ; N ¼ ; i;j i;j kþ1 kþ1 4s1 þa2 a2 þ4s1
kþ1=2 ¼ Fi;j
! kþ1=2 kþ1=2 kþ1=2 kþ1=2 kþ1=2 þ ci1;j Dkþ1=2 2 ciþ1;j 2ci;j 1 ciþ1;j ci1;j hz þ $ Dkþ1 h2r 2hr ri h i h i kþ1=2 1 kþ1 kþ1=2 kþ1=2 2 1 þ kDkþ1 1 þ 0:5sa2 Vi;j þ hz Dkþ1 Ni;j þ a3 gkþ1 ; ci;j
at every fixed i ¼ 2, 3,., N1 1 obtain kþ1=2 kþ1 kþ1 kþ1 ckþ1 ¼ 0; k 0: i;j1 Gi;j ci;j þ ci;jþ1 þ Fi;j
(19)
Find the concentration in grid points on the (k þ 1)th time layer in the form (20)
The sweep coefficients are: j ¼ 2; 3; .; N2 ;
ð15Þ
Find wkþ1 1;j from (8) (i ¼ 1, s ¼ 1/2) and substitute it into (15). In notations h i kþ1=2 kþ1=2 a1 1W1;j kþ1 kþ1=2 2 1 kþ1 1 ~ ; k ¼ 2hz s ; G1;j ¼ 2þkDkþ1 1þ V1;j ; V1;j ¼ akþ1 þ4s1 2
~ kþ1 kþ1=2 kþ1=2 1 W 2akþ1 1 1;j 4s1 a2 akþ1 2 w1;j kþ1 kþ1=2 ~ ¼ V ; N ¼ ; 1;j 1;j þ4s1 4s1 þakþ1 akþ1 2 2
2Dkþ1=2 2 kþ1=2 kþ1=2 kþ1=2 kþ1=2 F1;j ¼ ðhz =hr Þ c2;j 2c1;j þc0;j Dkþ1 h i h i kþ1=2 kþ1 kþ1=2 kþ1=2 c1;j þh2z D1 þkD1 þa3 gkþ1 ; kþ1 1þ0:5sa2 V1;j kþ1 N1;j we get
1 1 kþ1=2 kþ1 kþ1 ¼ Gkþ1 ; bkþ1 ¼ bkþ1 : Gkþ1 akþ1 i;j i;j1 þ Fi;j1 i;j i;j1 ai;j1 i;j1 ai;j1 Find the first coefficients from (19) at j ¼ 1 and from the condition cz jH=2 ¼ 0 : 1 ~ kþ1 ; bkþ1 ¼ Fkþ1=2 =2: 1þV akþ1 i;1 1;1 i;1 ¼ 1 kð2Dkþ1 Þ i;1 The immediate task is to find the values of ckþ1 i;N2 , i ¼ 1,., N1 1 which is necessary to launch the sweep method. Write down the approximation of the boundary condition (4) ðz ¼ HÞ : k 0; . h
2
kþ1=2 kþ1=2 ð0:5sÞ ¼ 0:5 bkþ1=2 q2i qkþ1 2i q2i i 2 Dkþ1=2 cz tkþ1=2 ; ri ; H bkþ1 qkþ1 Dkþ1 cz ðtkþ1 ; ri ; HÞ : 2i
ð21Þ
In the boundary node with a precision of up to Oðh2z Þ we have
kþ1=2 kþ1 kþ1 kþ1 ckþ1 ¼ 0; k 0: 1;j1 G1;j c1;j þ c1;jþ1 þ F1;j
(16) th
Find approximation of concentration on the (k þ 1) layer in the form
time
(17)
The sweep coefficients are: j ¼ 2; 3; .; N2 ; akþ1 1;j
1 ~ kþ1 ; bkþ1 ¼ Fkþ1=2 =2: 1 þ V ¼ 1 kð2D Þ akþ1 kþ1 1;1 1;1 1;1 1;1
kþ1 kþ1 ckþ1 ¼ akþ1 i;j i;jþ1 ci;jþ1 þ bi;jþ1 : j ¼ 0; 1; .; N2 1; k 0:
~ kþ1 ckþ1 1W akþ1 1;j 1 1;j
kþ1 kþ1 kþ1 ckþ1 1;j ¼ a1;jþ1 c1;jþ1 þ b1;jþ1 ; j ¼ 0; 1; .; N2 1; k 0:
1243
1 1
kþ1=2 kþ1 kþ1 kþ1 Gkþ1 ¼ Gkþ1 ; bkþ1 : 1;j ¼ b1;j1 þ F1;j1 1;j1 a1;j1 1;j1 a1;j1
Find the first coefficients from (16) at j ¼ 1 and from the condition cz jH=2 ¼ 0 :
m m 2hz cz ðtm ; ri ; HÞzcm i;N2 2 4ci;N2 1 þ 3ci;N2 ; m ¼ k þ 1=2; k þ 1:
(22)
Concentration values on the (k þ 1/2)th layer are already known. For the (k þ 1)th layer, substituting values ckþ1 i;N2 2 and from expression (17) (at i ¼ 1) and (20) (at i ¼ 2,., N1 1), ckþ1 i;N2 1 we get the approximation
1 h kþ1 kþ1 3 þ akþ1 cz ðtkþ1 ; ri ; HÞz ckþ1 i;N2 ai;N2 1 4 i;N2 þ bi;N2 1 2hz
i kþ1 þ akþ1 i;N2 1 4 bi;N2 : Write down the approximation in a more compact form:
1244
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 6 ( 2 0 1 1 ) 1 2 3 9 e1 2 4 7
Fig. 2 e The impact of the parameter b0 .
Fig. 1 e Diffusion flux approximations of IðtÞ.
kþ1 kþ1 kþ1 kþ1 Ah3 þ akþ1 i;N2 ai;N2 1 4 ; Bhbi;N2 1 þ ai;N2 1 4 bi;N2 ; 2hz cz ðtkþ1 ; ri ; HÞzAckþ1 i;N2 þ B:
ð23Þ
Expressions (22), (23) substitute into (21), denote skþ1 i;N2 ¼ y: bkþ1 2 4 Dkþ1 A y þ G ¼ 0; y þ þ sgkþ1 2hz g2kþ1 ! bkþ1=2 kþ1=2 2 3Dkþ1=2 4 c ckþ1=2 Gh 2 þ 2hz sgkþ1=2 i;N2 gkþ1=2 i;N2
1 kþ1=2 kþ1=2 þ Dkþ1 B þ Dkþ1=2 ci;N2 2 4ci;N2 1 : 2hz
1 m qm 2i ¼ gm ci;N2 ; m ¼ k þ 1=2; k þ 1;
Taking the physical meaning of the process into account we choose the positive root of the square equation with respect to y. The roots have different signs, at least at relatively small s (denote s 1). The approximation order of the boundary condition is Oðh2z þ s2 Þ. It agrees with the approximation order of the scheme in the bulk. Now find all values of ckþ1 i;j at i ¼ 0 and i ¼ N1, j ¼ 0, 1,., N2. Using boundary condition on the axis of the cylinder kþ1 kþ1 kþ1 ðcr jþ0 ¼ 0Þ, we get ckþ1 0;j ¼ ð4c1;j c2;j Þ=3. Values of cN1 ;j > 0 are defined from the equation
h i n kþ1=2 kþ1=2 kþ1=2 kþ1=2 ci;j wkþ1 ¼ 4s1 wi;j þ a1kþ1=2 1 Wi;j a2kþ1=2 wi;j i;j ~ kþ1 ckþ1 akþ1 þ 4s1 1 : 1W þ akþ1 1 2 i;j i;j ~ s (s ¼ k þ 1/2; k þ 1) is the following: we The correction of W i;j s ~ ¼ ws =wmax and proceed computing using the assign W i;j i;j scheme of the method of alternating directions until stabili~ s zWs (usually twoethree iterations). zation W i;j i;j
3.4.
Classical diffusion problem
For simplified (rougher) calculations one can suppose that diffusion is the only limiting process (physicalechemical processes on the surface are ignored); then the Dirichlet boundary-value problem with zero boundary concentrations is considered:
v2 c 1 vc v2 c vc ; r˛ð0;LÞ; z˛ð0;HÞ; t˛ð0; t Þ; þ $ þ ¼ DðtÞ vt vr2 r vr vz2 cr ðt; þ0; zÞ ¼ 0; cðt; L;zÞ ¼ 0; cðt; r;0Þ ¼ cðt; r; HÞ ¼ 0; cð0;r; zÞ ¼ c: Partial sum of series
bkþ1 kþ1 2 2 3Dkþ1 kþ1 Dkþ1 kþ1 cN1 ;j þ cN1 2;j 4ckþ1 cN1 ;j þ þ N1 1;j 2 2hr 2hr sgkþ1 gkþ1 2 kþ1=2 c ¼ 0; s gkþ1=2 N1 ;j which approximates the condition (3) at r ¼ L. Using ckþ1 i;j compute
Table 1 e Parameter values common for all curves.a b0 ¼ 1013 D0 ¼ 106 g0 ¼ 103
Eb ¼ 12 104 ED ¼ 5 104 Eg ¼ 0
T0 ¼ 293 Tcrit ¼ 700 Te ¼ 773
a1 ¼ 0:1 a2 ¼ 0:2 a3 ¼ 103
L ¼ 4 103 H ¼ 2 102 c ¼ 1023
a ½b0 ; D0 ¼ m2 s1 , ½Eb ; ED ; Eg ¼ J mol1 , ½g0 ¼ m1 , ½T0 ; Tcrit ; Te ¼ K, ½a1 ; a2 ; a3 ¼ s1 , ½L; H ¼ m, ½c ¼ m3 .
Fig. 3 e The impact of the parameter Eb .
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 6 ( 2 0 1 1 ) 1 2 3 9 e1 2 4 7
Fig. 4 e The impact of the parameter D0 .
8 9 Zt < = m0
np cðt; r; zÞ ¼ Anm exp lnm DðsÞds $J0 m r $sin z ; : ; L H n¼1 m¼1 0
np2 m0 2 4 1 þ ð1Þnþ1 c ; lnm h þ m Anm h L H pnm0m J1 m0m N X N X
is used as approximation of the solution. This sum is a generalized solution since boundary conditions are inconsistent. Here, DðtÞhDðTðtÞÞ, J0, J1 are zero-order and first-order Bessel functions of the first kind, respectively, m0m are sequential zeros of the function J0(m). H flux through the surface is: Z Dcn ds ¼ c
IðtÞ ¼
i h nþ1 N X N 8HDðtÞ 1 þ ð1Þ X
np n¼1 m¼1 " 2 # 1 þ ð1Þnþ1 pL $expf.g: þ ð1Þnþ1 2n n Hm0m
1245
Fig. 6 e The impact of the parameter g0 .
time is not so important: interesting effects take place in the neighborhood of the TDS-peak “top”. How many terms of a series should we include in the calculations for appropriate curve approximation (5%)? The series converges slowly. For instance, at D0 ¼ 2 103 m2 s1, ED ¼ 6 104 J mol1 (Dz2 107 , T ¼ 773 K ), c ¼ 1023 m3 one should take 144 items of a series (n ¼ 1, 3, 5, 7, 9, 11, 13, 15; 1 m 18). The series equals zero when n is even. Fig. 1 shows the approximation of I(t) with partial sums (n, m) (N, M ) ¼ (15, 18), (5 10), (3, 3), (1, 1) with respect to the decrease of the maximum. It is insufficient to include several items, the peak is unique (in the concrete situation it is necessary to take several comparable diffusion channels with different D( j ) into consideration).
S
4. Algorithm and results of numerical modelling
One often takes only the first item (n ¼ m ¼ 1) into account, believing that the rest of the items are quickly damped. Yet, in the problem under consideration the asymptotic form versus
Describe the computing algorithm step-wise. We define the values of L, H, D0, ED, b0, Eb, g0, Eg, c, wmax, g, T0, Te, Tcrit, a1, a2,
Fig. 5 e The impact of the parameter ED .
Fig. 7 e The impact of the parameter Eg .
1246
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 6 ( 2 0 1 1 ) 1 2 3 9 e1 2 4 7
Fig. 9 e The impact of the defects of two types. Fig. 8 e Comparison of JðtÞ and IðtÞ.
a3, s, r, ~c. Transition from the kth layer to the (k þ 1)th layer is realized in two steps. Step I: the computing algorithm on the (k þ 1/2)th layer ðk 0Þ is the following. ~ from equation (9) at 1. Compute the values of W i;j s ¼ k þ 1=2. 2. According to expressions (10), (11) compute the sets of kþ1=2 kþ1=2 , bi;j . coefficients ai;j 3. Determine the values of concentration in boundary nodes, kþ1=2 solving the square equation with respect to y ¼ sN1 ;j > 0. 4. Using formula (11) find the approximate values of concentration in all internal nodes. kþ1=2 in 5. Using boundary conditions determine the values of ci;j boundary nodes at j ¼ 0 and j ¼ N2 , i ¼ 0; 1; .; N1 . 6. From the second equation of the difference scheme in the bulk compute values of concentration in defects of the reversible capture type:
6. From the second equation of the difference scheme in the bulk obtain: h i N kþ1=2 kþ1=2 kþ1=2 kþ1 wi;j ¼ kþ1 ; N ¼ 4s1 wi;j þ a1kþ1=2 1 Wi;j ci;j a2 þ 4s1 kþ1=2 ~ kþ1 ckþ1 : 1W a2kþ1=2 wi;j þ akþ1 i;j 1 i;j
kþ1=2
kþ1=2 ¼ wi;j
h i kþ1=2 k ~ kþ1=2 ckþ1=2 1 W 4s1 wki;j þak1 1Wi;j cki;j ak2 wki;j þa1 i;j i;j
7. Correction
a2kþ1=2 þ4s1 of
~ kþ1=2 W i;j
is
the
following:
one
:
assign
~ kþ1=2 ¼ wkþ1=2 =wmax and proceed as in items 2e6 until W i;j i;j ~ kþ1=2 zWkþ1=2 (twoethree iterations). stabilization W i;j i;j Step II: the computing algorithm on the (k þ 1)th time layer ðk 0Þ is the following. ~ from equation (9) at s ¼ k þ 1. 1. Compute the values of W i;j 2. According to expressions (16), (17), (19), (20) compute coefkþ1
kþ1 kþ1 , bkþ1 ficients akþ1 1;j , ai;j , bi;j . 1;j
3. Determine the values of concentration in boundary nodes, solving the square equation with respect to y ¼ ckþ1 i;N2 > 0. 4. Using formulas (17), (20) find the approximate values of concentration in internal nodes. at i ¼ 0 and 5. From boundary conditions determine ckþ1 i;j i ¼ N1 , j ¼ 0; 1; .; N2 .
7. Correction
of
~ kþ1 W i;j
is
the
following:
one
assign
~ kþ1 ¼ wkþ1 =wmax and proceed as in items 2e6 until stabiW i;j i;j kþ1 ~ lization W i;j zWi;j . kþ1
4.1.
Results of numerical modelling
Desorption flux curves JðtÞ are represented in the figures. Varied coefficient values correspond to the curves in the following way: the first value is for the curve with the highest maximum, the last one is for that with the lowest maximum. The small circle shows the time instant when the initial surface hydrogen is exhausted. The algorithm developed for numerical modelling of hydrogen thermodesorption is applicable in a wide range of parameters. In our figures however, we limit ourselves to concrete values of the parameters presented below. We define g ¼ 1023 m3 , wmax ¼ 1023 m3 . The default parameters, which are common for all curves, are taken from the Table 1. Fig. 2 illustrates the influence of the coefficient b0 on the desorption flux curve. Changes are: b0 ¼ 1012 ; 1013 ; 1014 m2 s1 . The influence of the activation energy Eb can be seen in Fig. 3, changes in order of decrease of the maximum are: Eb 104 ¼ 9; 10; 12 J mol1 . Changes for Fig. 4 are: D0 ¼ 5 105 ; 105 ; 106 m2 s1 . Fig. 5 shows the influence of the parameter ED : ED 104 ¼ 4; 5; 6. The variations of g0 ¼ 50; 100; 200 m1 are seen in Fig. 6. Fig. 7 shows the dependence of the desorption curve on the value of Eg 103 ¼ 1; 0; 1. Fig. 8 shows two flux curves: the left-hand curve is the desorption flux for the accepted boundary-value problem with dynamic boundary conditions, the right-hand curve is the diffusion flux for the degenerate Dirichlet problem. Parameter values are: D0 ¼ 2 103 , ED ¼ 6 104 , b0 ¼ 105 , Eb ¼ 12 104 , g0 ¼ 103 , Eg ¼ 0, a10 ¼ 0, a3 h0, g ¼ 3 1023 m3 ,
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 6 ( 2 0 1 1 ) 1 2 3 9 e1 2 4 7
wmax ¼ 1024 m3 . The flux curve taking the defects of two types into account is illustrated in Fig. 9 (reversible capture and decomposition).
Acknowledgements The work was supported by the Russian Foundation for Basic Research (grant 09-01-00439) and by the Russian Science Support Foundation.
[4]
[5]
[6]
references
[1] Kolachev BA. Hydrogen embrittlement of metals. Moscow: Metallurgy; 1985 [in Russian]. [2] Kunin LL, Golovin AI, Surovoi UI, Hohrin VM. Problems of metals degassing. Moscow: Nauka; 1972 [in Russian]. [3] Polyanskiy AM, Polyanskiy VA, Popov-Diumin DB, Kozlov EA. New measuring complex for absolute determination of
[7]
[8]
1247
hydrogen content in hydrogen energetics materials. Int J Altern Energy Ecology 2006;6:29e31. Gabis IE, Kompaniets TN, Kurdyumov AA. In: Zakharov AP, editor. Interactions of hydrogen with metals. Moscow: Nauka; 1987. p. 177e206 [in Russian]. Polyanskiy AM, Polyanskiy VA, Yakovlev Yu A. Methodic of investigation of the binding energies of hydrogen in solid state on the basis of hydrogen analyzer AV-1. Proc. III Int. Conf. “Interaction of hydrogen isotopes with structural materials”. Sarov, RFNCeVNIIEF; 2008. p. 346e53. Zaika YuV, Rodchenkova NI. Modelling of diffusion TDSspectrum peak of dehydriding with size reduction and heat absorption effects. In: Baranowski B, et al., editors. Carbon nanomaterials in clean energy hydrogen systems, NATO science for peace and security series e C: environmental security. Springer; 2008. p. 863e78. Zaika YuV, Rodchenkova NI. Boundary-value problem with moving bounds and dynamic boundary conditions: diffusion peak of TDS-spectrum of dehydriding. Appl Math Model 2009; 33:3776e91. Samarskij AA. Introduction to difference schemes theory. Moscow: Nauka; 1971 [in Russian].
Available at www.sciencedirect.com
journal homepage: www.elsevier.com/locate/he
Numerical modelling of hydrogen desorption from cylindrical surface Natalia I. Rodchenkova*, Yury V. Zaika Institute of Applied Mathematical Research, Karelian Research Centre of the Russian Academy of Sciences, Pushkinskaya 11, 185910 Petrozavodsk, Russian Federation
article info
abstract
Article history:
We consider degassing of cylindrical sample containing dissolved hydrogen. The experi-
Received 25 June 2010
ment is made by the thermodesorption method. In the corresponding boundary-value
Accepted 25 June 2010
problem with nonlinear dynamic boundary conditions physicalechemical processes in the
Available online 17 September 2010
bulk and on the metal surface are taken into account: diffusion, desorption, solution, and capture by defects. Computational algorithm for desorption flux modelling is developed on
Keywords: Hydrogen thermodesorption
the basis of difference approximations. The results of numerical modelling are presented. ª 2010 Professor T. Nejat Veziroglu. Published by Elsevier Ltd. All rights reserved.
Mathematical modelling
1.
Introduction
A high concentration of hydrogen in metal leads to hydrogen embrittlement [1,2]. Common metallurgical concentrations of dissolved hydrogen range from 0.1 to 100 ppm. The authors of paper [3] have developed the hydrogen analyzer AV-1. The device allows to measure the concentration of hydrogen atoms (H) in a solid sample using industrial laboratory facilities. A cylindrical sample is placed into the quartz vacuum extractor. The extractor is exposed to a given temperature in the furnace. The contact between the sample and the extractor walls is point-like, the thermal conductivity of quartz is negligible, wherefore the heat transfer is determined by radiation. As the sample is heated, atomic hydrogen diffuses inside the bulk and is desorbed from the surface in the molecular form. The extraction curve (measured by a mass-spectrometric hydrogen analyzer) is recorded. The curve is used for further processing (in particular, the kinetic parameters of models are estimated). The dependence of the desorption flux on the temperature at uniform heating (the TDS-spectrum) usually has several peaks. According to [4], in addition to the diffusion, other surface processes may affect
the desorption of hydrogen (e. g. the capture of hydrogen atoms by cracks, microcavities; inclusion of hydrogen into hydride phases). The experimental approaches to the problems of hydrogen materials science, especially in the case of heavy isotopes deuterium and tritium, are rather complicated and expensive. Therefore, the role of mathematical modelling in hydrogen materials science is very important. The paper is devoted to the development of mathematical tools for reliable experimental measurements of the hydrogen content in various materials. In this work, we will use parameters suitable for aluminium and some of its alloys when numerically modelling the desorption flux.
2.
Mathematical model
2.1.
Heating equation
The sample has the form of a cylinder with typical dimensions (GOST 21132.1-98): radius of the base is L ¼ 4 103 m, height is H ¼ 2102 m. Concentration of diffusing atomic hydrogen at the initial time instant cð0; r; zÞ ¼ c is constant (it forms during
* Corresponding author. Tel.: þ7 8142766312; fax: þ7 8142766313. E-mail addresses: [email protected] (N.I. Rodchenkova), [email protected] (Yu.V. Zaika). 0360-3199/$ e see front matter ª 2010 Professor T. Nejat Veziroglu. Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.ijhydene.2010.06.121
1240
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 6 ( 2 0 1 1 ) 1 2 3 9 e1 2 4 7
Nomenclature L, H S, V t ts T Te tcrit, Tcrit r s ~c l g
radius, height of the cylinder (m) surface area and capacity (m2, m3) time (s) time instant when surface H is exhausted (s) temperature (K) temperature of extractor wall (K) critical time and temperature (s, K) bulk density of the metal (kg/m3) Stefan-Boltzmann constant (J/s/m2/K4) specific heat capacity (J/kg/K) thermal conductivity (J/s/m/K) quick dissolution coefficient (1/m)
the material production). When necessary, it is possible to take account of the decrease in H concentration in the near-surface layer (for instance, as the result of mechanical and thermal pretreatment) without fundamental changes in the numerical algorithm. If the sample heating is uniform (sufficiently slow, T ¼ T(t), ½T ¼ K), the temperature dynamics can be described by the differential equation of the heat balance [5]: dT sS 7 105 ðT þ 64:3Þ T4e T4 ; ¼ dt ~crV
T0 ¼ Tð0Þ ¼ 293 K:
Here, Te ¼ const is the temperature of the extractor wall, S ¼ 2pLH, V ¼ pL2H are the surface area and the cylinder capacity, r ¼ 2.71 103 kg m3 is the bulk density, s ¼ 5.67 108 J s1 m2 K4 is the Stefan-Boltzmann constant, ~c ¼ 1:15 103 J kg1 K1 is the specific heat capacity. Since the assumption about uniform heating does not have adequate accuracy for all values of Te, we consider an alternative distributed model. We suppose that in the case of tubular geometry of the extractor heating mainly proceeds through the lateral surface. Thus, we speak of a lower estimate of the heating of the centre of the sample. The considered equation from [5] is majorized from above. The radiosymmetrical model is: T(0, r) ¼ T0 < Te, r˛½0; L,
J hydrogen desorption flux (1/s) k thermal diffusivity (m2/s) heat absorption coefficient 3t a emissivity (W/m2/K4) R universal gas constant (J/mol/K) r, z cylindrical coordinates (m) c H concentration in metal (1/m3) w H concentration in defects (1/m3) g H concentration in hydride phase (1/m3) c; w; g initial concentrations (1/m3) surface concentrations (1/m2) qi defect coefficients (1/s) ai E activation energy (J/mol) D, b diffusion and desorption (m2/s)
2.2.
Diffusion model with defects
Consider the boundary-value problem of thermodesorption for the cylinder taking into consideration diffusion in the bulk, H capture by defects, penetration from the bulk to the surface and desorption: 2 vc v c 1 vc v2 c þ $ þ 2 a1 ½1 Wcðt; r; zÞ ¼DðTÞ 2 vt vr r vr vz (1) þ a2 wðt; r; zÞ þ a3 gðtÞ; vw w ; ¼ a1 ðTÞ½1 Wc a2 ðTÞw; Wh vt wmax (2) r˛ð0; LÞ; z˛ð0; HÞ; t˛ð0; t Þ;
2 vT v T 1 vT vT vT þ ¼ 0; l ¼ a T4e T4 ðt; LÞ : ¼k $ ; vt vr2 r vr vr r¼þ0 vr L
dg ¼a3 gðtÞ; gðtcrit Þ¼g0gðtÞ¼gexpfðtcrit tÞa3 g; a3 ¼0; T
Here, l is the heat conduction coefficient (for aluminium in the range of T˛½300; 800 lz236 J s1 m1 K1 ); k ¼ l ð~crÞ1 is the thermal diffusivity; the finishing time of computing t* is determined by the stationary distribution Tðt; 0ÞzTe , t > t*; a ¼ s3t, 3t ¼ 7 105(T þ 64.3) is the heat absorption coefficient. The auxiliary problem of numerical modelling of heating is the following: it is necessary to quickly estimate how much the distribution T(t, r) differs from the uniform heating T(t) under given Te, L, H, and thermal and physical characteristics of the material. For instance, when we use T0 ¼ 293 K, Te ¼ 773 K and the L, H stated above the assumption about uniform heating of aluminium sample is true: the difference T(t) T(t, 0) is within a hundredth of a degree. We obtain this result in the distributed model where heat absorption by butt ends is disregarded. Peak temperature is reached in 2.2 h, heating being practically uniform within an hour. Hereinafter, suppose that the sample heating is uniform.
Here, c(t, r, z) is the concentration of corpuscular (diffusion movable) hydrogen in the metal; w(t, r, z) is H concentration in the defects of the material physicalechemical structure, the capture is reversible (for instance, microcavities); wmax is the peak concentration of reversible capture; g(t) is H concentration in traps, which begin to release the hydrogen only when some critical temperature Tcrit is reached (it is typical for the inclusion of hydride phases); ai are the coefficients of absorption and desorption of H by traps (a3 > 0 only at T Tcrit); q1(t, z), q2(t, r), q3(t, r) are the surface concentrations (on the lateral surface of the cylinder and on butt ends); c is the initial concentration c(0,r,z) in the solid sample; g is the fit coefficient for concentrations of hydrogen atoms in the bulk and on the surface (quick dissolution coefficient); D, b are the diffusion and the desorption coefficients. Suppose that ai > 0 are constant within the considered temperature range ðT˛½300; 800Þ. Changes in the
1241
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 6 ( 2 0 1 1 ) 1 2 3 9 e1 2 4 7
time-dependent case of ai ðtÞhai ðTðtÞÞ are insignificant. For practical purposes we consider capture in a simple integrated form. Specifications of defects geometry and their distributions in the sample essentially complicate the model. To simplify the computations we consider one generalized trap with reversible capture and one generalized trap with decomposition. It is possible to consider several ðjÞ ðkÞ traps with their individual ai and Tcrit . We assume that diffusion and desorption coefficients depend on the temperature in an Arrhenius way: D(T ) ¼ D0exp{ED/[RT]}, b(T ) ¼ b0exp{Eb/[RT]}, ED,b are activation energies. We use a shorthand notation D(t) ¼ D(T(t)), b(t) ¼ b(T(t)). The time segment [0, t*] is determined by total degassing of the sample: JðtÞz0, t t*. The conditions cr(t, þ0, z) ¼ 0, cz(t, r, H/ 2) ¼ 0 are consequences of symmetry. A more precise model of dissolution on the surface (considering the lateral surface for clarity) has a form of fluxes balance: 1 kþ ðTÞcðt; L; zÞ 1 q1 ðt; zÞq1 max k ðTÞq1 ðt; zÞ 1 cðt; L; zÞcmax ¼ DðTÞcr ðt; L; zÞ: However, when diffusion is considerably slower than dissolution (the temperature is not very low) and concentrations are low, we obtain the condition of quick dissolution czgq, where g ¼ k/kþ. If the surface is isotropic (in the sense Ek zEkþ ), then the parameter g weakly depends on T. Below we use the notation Eg conditionally: it is not the activation energy (the difference Ek Ekþ may be negative). For traps which begin to release hydrogen only when certain critical temperature is reached (traps of the hydride phase type) we did the following: we took into account only the trap capacity and decomposition rate. Modelling of the dehydriding process is an independent intricate problem, which results in nonlinear boundary-value problems with free phase boundaries and with conditions of the Stefan type. The variant with bulk desorption is considered in Ref. [6], the variant with surface desorption is presented in Ref. [7]. Since the initial data are symmetrical we have q3 ¼ q2, and so below we construct the difference approximation only for a half of the cylinder ðz˛½H=2; HÞ with corresponding boundary conditions ðcz jH=2 ¼ 0; q_ 2 ¼ .Þ. For the defect with reversible capture (microcavity) the constant w is determined according to (2) ðv=vt ¼ 0Þ : a1 ðT0 Þ½1 w=wmax c a2 ðT0 Þw ¼ 0. For the trap of the hydride phase inclusion type the values of g ¼ const, Tcrit, a3 are set using information about specific chemistry of the hydride. The presence of derivatives q_ i (accumulation on the surface) corresponds to the possibility of H migration along the surface until H atoms form H2 molecules which are desorbed from the surface. Thus, the model aimed to analyze the problem of the dynamics of low common concentrations of hydrogen (without preliminary saturation in laboratory). The problem has an applied context, so the authors apply minimal mathbased environment to describe basic processes. Further detailed elaboration leads to an increase in the number of parameters. The inverse problem of the parameters estimation becomes complicated.
2.3.
Specification of the problem statement
The aim of the paper is to develop the difference scheme and computational algorithm for modelling of hydrogen desorption flux from cylindrical sample: Z JðtÞ ¼ 4pbðtÞ L
Z
H H=2
q21 ðt; zÞdz
þ
L 0
! rq22 ðt; rÞdr
:
The H2 molecule is desorbed from the surface, but we compute the desorption flux in atoms ([J] ¼ 1/s). The criterion of computation correctness is the material balance: Z pL2 Hðc þ w þ gÞ þ 4p L Z ¼ 4p
Z
H
Z þ 4p L
Z q1 ð0; zÞdz þ
H=2 L
dz H=2
H
! rq2 ð0; rÞdr
0
r½cðt; r; zÞ þ wðt; r; zÞ þ gðtÞdr
0 H
L
Z q1 ðt; zÞdz þ
H=2
L
! rq2 ðt; rÞdr
0
Z þ
t
JðsÞds: 0
Under uniform heating, it is convenient to consider the TDS-spectrum (the curve J ¼ J(T )) together with the flux dependence on t. TDS-spectrum usually has several peaks. It is believed the first peak corresponds to surface hydrogen. But one should have caution: while the surface hydrogen is desorbed, hydrogen atoms penetrate from the bulk to the surface. The problem of estimating the corresponding correction is topical. Let ts be the time instant when the surface hydrogen is exhausted. At the numerical modelling stage the time instant ts is determined from Z
ts
Z JðsÞdsz4p L
0
H H=2
Z q1 ð0; zÞdz þ
L
! rq2 ð0; rÞdr :
0
Numerical modelling allows to choose a region on the extraction curve which corresponds to the initial quantity of the surface hydrogen, to evaluate the value of the activation energies of diffusion and desorption, to estimate the parameters of reversible capture and the parameters of hydride phases decomposition.
3. Difference approximation of boundaryvalue problem 3.1.
Difference scheme construction
Follow the standard method [8], consider a spatial grid Uh ¼ ri ; zj : ri ¼ ihr ; i ¼ 0; 1; .; N1 ¼ ½L=hr ; zj ¼ jhz ;
j ¼ 0; 1; .; N2 ¼ ½ðH=2Þ=hz and a time grid us¼{tk ¼ ks, k ¼ 0, 1,., K ¼ [t*/s]}. Denote approximate values of the bulk concentration c(tk, ri, zj) by cki;j . Similarly, wki;j zwðtk ; ri ; zj Þ; gk zgðtk Þ; aks ¼ as ðtk Þ, Dk ¼ D(tk), where ðri ; zj Þ˛Uh , tk ˛us . For equation (1) examine the implicit difference scheme of the alternating directions method, called wedge-reed (PeacemaneRachford) scheme, while for equation (2) consider the scheme with weights. Transition from the kth layer to the (k þ 1)th layer is realized in two steps. At the first kþ1=2 are determined from the step intermediate values of ci;j equations set
1242
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 6 ( 2 0 1 1 ) 1 2 3 9 e1 2 4 7
kþ1=2 ci;j cki;j
0:5s
kþ1=2 kþ1=2 kþ1=2 ciþ1;j 2ci;j þ ci1;j
kþ1=2 kþ1=2 1 ciþ1;j ci1;j ¼ Dkþ1=2 þ $ h2r 2hr ri k k k ci;jþ1 2ci;j þ ci;j1 kþ1=2 kþ1=2 ~ þ Dk ci;j a1kþ1=2 1 W i;j h2z
!
kþ1=2 þ a2kþ1=2 wi;j þ a3 gkþ1=2 ;
kþ1=2 wi;j wki;j
0:5s
ð5Þ
h i k ¼ð1 sÞ ak1 1 Wi;j cki;j ak2 wki;j ~ kþ1=2 ckþ1=2 akþ1=2 wkþ1=2 : þ s akþ1=2 1 W 1
kþ1=2 kþ1=2 kþ1=2 kþ1=2 Ai ci1;j Gi;j ci;j þ Bi ciþ1;j þ Fki;j ¼ 0; k 0:
i;j
i;j
2
i;j
kþ1=2 ci;j
0:5s
kþ1=2 ciþ1;j
þ
kþ1=2 ci1;j
kþ1=2 ciþ1;j
kþ1=2 ci1;j
1 þ $ 2hr ri kþ1 kþ1 ckþ1 i;jþ1 2ci;j þ ci;j1 kþ1 ~ kþ1 ckþ1 1 W þ Dkþ1 a i;j 1 i;j h2z
¼Dkþ1=2
kþ1=2 2ci;j h2r
kþ1 kþ1 þ akþ1 ; 2 wi;j þ a3 g
kþ1=2
wkþ1 i;j wi;j 0:5s
ð6Þ
!
i;j
At r / þ0 we have cr =r ¼ ðcr ðt; r; zÞ cr ðt; 0; zÞÞ=rzcrr . Find the first coefficients from approximation of the equation
1 kþ1=2 ~ kþ1=2 ; bkþ1 ¼ Fk =4: 1þV cr jþ0 ¼ 0 : a1;j ¼ 1 k 4Dkþ1=2 1;j 1;j 1;j which is The immediate task is to find a value of cNkþ1=2 1 ;j necessary to launch the sweep method. Write down the approximation of the boundary condition (3) (r ¼ L) : (k 0),
. h
2 ð0:5sÞ ¼ 0:5 bk qk1j
2 i kþ1=2 Dk cr tk ; L; zj bkþ1=2 q1j Dkþ1=2 cr tkþ1=2 ; L; zj :
kþ1=2
q1j
qk1j
(9)
Sweep method with respect to the radius r
ð12Þ
In the boundary node with a precision of up to Oðh2r Þ we have m m 2hr cr tm ; L; zj zcm N1 2;j 4cN1 1;j þ 3cN1 ;j ; m ¼ k; k þ 1=2:
In the standard notations y_ ¼ f ðt; yÞ, ys ¼ ys1/2 þ {fs1/2 þ fs}s/4, the equation (9) is the symmetrical scheme for equation (2) normalized to wmax with the time-fixed function c ¼ c(ts1/2, r, z). We explain below the iterative procedure of more precise ~ s . For distinctness define s ¼ 1/2. The approxdefinition of W i;j imation order is Oðs2 þ h2r þ h2z Þ [8].
(13)
Concentration values on the kth layer are already known. For and cNkþ1=2 from the (k þ 1/2)th layer, substituting values cNkþ1=2 1 2;j 1 1;j expression (11), we have
1 h 3 þ aNkþ1=2 aNkþ1=2 cr tkþ1=2 ; L; zj z cNkþ1=2 þ bNkþ1=2 ;j 1;j 4 1 1 1 ;j 1 1;j 2hr i
: þ aNkþ1=2 4 bNkþ1=2 1 1;j 1 ;j Write down the approximation in a more compact way:
kþ1=2 kþ1=2 kþ1=2 kþ1=2 kþ1=2 Ah3 þ aN1 ;j aN1 1;j 4 ; BhbN1 1;j þ aN1 1;j 4 bN1 ;j ; þB: 2hr cr tkþ1=2 ; L; zj zAcNkþ1=2 1 ;j
ð14Þ
th
Consider a transition from the k layer to the (k þ 1/2) layer. kþ1=2 from the equation (6) and substitute it into (5). Express wi;j In notations Ai ¼ 1hr ð2ri Þ1 ; Bi ¼ 1þhr ð2ri Þ1 ; k ¼ 2h2r s1 ; h i k ak1 1Wi;j kþ1=2 kþ1=2 k ~ ; Vi;j ¼ 2þkD1 ¼ kþ1=2 ; Gi;j kþ1=2 1þ Vi;j þ4s1 a2 kþ1=2 1 ~ kþ1=2 1 W 2a1 i;j 4s ak2 wki;j kþ1=2 k ~ V ¼ ; N ¼ ; i;j i;j 4s1 þa2kþ1=2 a2kþ1=2 þ4s1
Dk 2 Fki;j ¼ ðhr =hz Þ cki;jþ1 2cki;j þcki;j1 Dkþ1=2 h i h i kþ1=2 k 1 k cki;j þh2r D1 þkDkþ1=2 1þ0:5sa2kþ1=2 Vi;j Ni;j þa3 gkþ1=2 ; kþ1=2 a2 at every fixed j ¼ 1, 2,., N2 1 obtain
. . kþ1=2 kþ1=2 kþ1=2 kþ1=2 ¼Bi1 Gi1;j Ai1 ai1;j ; bi;j ¼ Ai1 bi1;j þ Fki1;j
kþ1=2 kþ1=2 ; i ¼ 2; 3; .; N1 : Gi1;j Ai1 ai1;j
on the (k þ 1/2)th time layer for i ¼ 1 and using the condition
h i
kþ1=2 kþ1=2 kþ1=2 ci;j ¼ð1 sÞ a1kþ1=2 1 Wi;j a2kþ1=2 wi;j ~ kþ1 ckþ1 akþ1 wkþ1 : 1 W ð8Þ þ s akþ1 i;j 1 2 i;j i;j
th
kþ1=2
ai;j
ct ¼ Dð2crr þ czz Þ a1 ½1 Wc þ a2 w þ a3 g
Here, in order to have a possibility to use the sweep method on s the sth time layer (s ¼ k þ 1/2; k þ 1), the unknown Wi;j is replaced by its approximation from the linear equation with ~s respect to W
3.2.
(11)
Sweep coefficients are:
ð7Þ
h i n ~ s ¼ Ws1=2 þ s as1=2 1 Ws1=2 cs1=2 w1 as1=2 Ws1=2 W i;j 2 max i;j i;j i;j i;j 4 1 h i o s s ~ : ~ cs1=2 w1 as W þ as1 1 W i;j i;j i;j max 2
Values at the initial time instant (on zero layer) are known: c0i;j ¼ c ¼ const. According to the sweep method, search approximate values of the concentration in grid points on the (kþ1/2)th time layer in the form kþ1=2 kþ1=2 kþ1=2 kþ1=2 ci;j ¼ aiþ1;j ciþ1;j þ biþ1;j ; i ¼ 0; 1; .; N1 1; k 0:
kþ1=2 At the second step, using the resultant ci;j , find ckþ1 from i;j the system
ckþ1 i;j
(10)
Expressions (13), (14) substitute into (12), denote sNkþ1=2 ¼ y: 1 ;j 1 m qm 1j ¼ gm cN1 ;j ; m ¼ k; k þ 1=2; # " 2 bkþ1=2 2 Dkþ1=2 A 4 bk y þ G ¼ 0; Gh 2 ckN1 ;j y þ þ 2hr sgkþ1=2 g2kþ1=2 gk
3Dk 4 1 Dkþ1=2 B þ Dk ckN1 2;j 4ckN1 1;j : ckN1 ;j þ þ 2hr sgk 2hr
This is a square equation with respect to variable y. The >0 roots have different signs (s 1). We choose y ¼ sNkþ1=2 1 ;j taking into account the physical meaning of the process. Approximation order of the boundary condition is Oðh2r þ s2 Þ. It agrees with the approximation order of the scheme in the bulk.
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 6 ( 2 0 1 1 ) 1 2 3 9 e1 2 4 7
kþ1=2
Now find all values of ci;j at j ¼ 0 and j ¼ N2, i ¼ 0, 1,., N1. Using the symmetry condition cz jH=2 ¼ 0, we have kþ1=2 kþ1=2 kþ1=2 kþ1=2 ¼ ð4ci;1 ci;2 Þ=3. Values ci;N2 > 0 (s, h are small) are ci;0 unambiguously defined from the square equation which approximates the boundary condition (4) at z ¼ H: # " bkþ1=2 kþ1=2 2 3Dkþ1=2 kþ1=2 2 ci;N2 þ þ c 2hz sgkþ1=2 g2kþ1=2 i;N2 2 Dkþ1=2 kþ1=2 kþ1=2 þ ci;N2 2 4ci;N cki;N2 ¼ 0: 2 1 2hz sgk kþ1=2 , compute the concenWhen we know all values of ci;j tration of w:
h i kþ1=2 k wi;j ¼ 4s1 wki;j þ ak1 1 Wi;j cki;j ak2 wki;j 1 ~ kþ1=2 ckþ1=2 akþ1=2 þ 4s1 þ a1kþ1=2 1 W : 2 i;j i;j
Sweep method with respect to the variable z
3.3.
Since there is a singularity in cylindrical coordinates at r / þ0, we make a transition from the (k þ 1/2)th layer to the (k þ 1)th layer in two steps. First step is: i ¼ 1, r / þ0 (approximate cr =rzcrr ), we realize the sweep method for the equation ct ¼ D(2crr þ czz) a1[1 W]c þ a2w þ a3g. Second step is: i ¼ 2,., N1 1, r > 0, realize the sweep method for the equation ct ¼ (crr þ cr/r þ czz) a1[1 W]c þ a2w þ a3g. The equation for w remains the same. Equation (7) for i ¼ 1 now is kþ1=2 ckþ1 1;j c1;j
0:5 s
¼2Dkþ1=2 þ þ
kþ1=2 kþ1=2 kþ1=2 c2;j 2c1;j þ c0;j
h2r
kþ1 2ckþ1 1;j þ c1;j1 Dkþ1 2 hz kþ1 akþ1 w þ a gkþ1 : 3 2 1;j
ckþ1 1;jþ1
Difference approximation of equation (7) for i ¼ 2,., N1 1 is: kþ1=2 ckþ1 i;j ci;j
kþ1=2 kþ1=2 1 ciþ1;j ci1;j þ $ 2 0:5s hr 2hr ri kþ1 kþ1 ckþ1 2c þ c i;jþ1 i;j i;j1 ~ kþ1 ckþ1 1W þ Dkþ1 akþ1 i;j 1 i;j h2z
¼ Dkþ1=2
kþ1=2 kþ1=2 kþ1=2 ciþ1;j 2ci;j þ ci1;j
kþ1 kþ1 þ akþ1 : 2 wi;j þ a3 g
!
ð18Þ
Find wkþ1 from (8) (s ¼ 1/2) and substitute it into (18). In i;j notations h i kþ1=2 a1kþ1=2 1Wi;j kþ1 kþ1=2 2 1 kþ1 1 ~ ; Vi;j ¼ ; k ¼ 2hz s ; Gi;j ¼ 2þkDkþ1 1þ V i;j akþ1 þ4s1 2
~ kþ1 kþ1=2 1 W 2akþ1 1 i;j 4s1 a2kþ1=2 akþ1 2 wi;j kþ1 kþ1=2 ~ ¼ V ; N ¼ ; i;j i;j kþ1 kþ1 4s1 þa2 a2 þ4s1
kþ1=2 ¼ Fi;j
! kþ1=2 kþ1=2 kþ1=2 kþ1=2 kþ1=2 þ ci1;j Dkþ1=2 2 ciþ1;j 2ci;j 1 ciþ1;j ci1;j hz þ $ Dkþ1 h2r 2hr ri h i h i kþ1=2 1 kþ1 kþ1=2 kþ1=2 2 1 þ kDkþ1 1 þ 0:5sa2 Vi;j þ hz Dkþ1 Ni;j þ a3 gkþ1 ; ci;j
at every fixed i ¼ 2, 3,., N1 1 obtain kþ1=2 kþ1 kþ1 kþ1 ckþ1 ¼ 0; k 0: i;j1 Gi;j ci;j þ ci;jþ1 þ Fi;j
(19)
Find the concentration in grid points on the (k þ 1)th time layer in the form (20)
The sweep coefficients are: j ¼ 2; 3; .; N2 ;
ð15Þ
Find wkþ1 1;j from (8) (i ¼ 1, s ¼ 1/2) and substitute it into (15). In notations h i kþ1=2 kþ1=2 a1 1W1;j kþ1 kþ1=2 2 1 kþ1 1 ~ ; k ¼ 2hz s ; G1;j ¼ 2þkDkþ1 1þ V1;j ; V1;j ¼ akþ1 þ4s1 2
~ kþ1 kþ1=2 kþ1=2 1 W 2akþ1 1 1;j 4s1 a2 akþ1 2 w1;j kþ1 kþ1=2 ~ ¼ V ; N ¼ ; 1;j 1;j þ4s1 4s1 þakþ1 akþ1 2 2
2Dkþ1=2 2 kþ1=2 kþ1=2 kþ1=2 kþ1=2 F1;j ¼ ðhz =hr Þ c2;j 2c1;j þc0;j Dkþ1 h i h i kþ1=2 kþ1 kþ1=2 kþ1=2 c1;j þh2z D1 þkD1 þa3 gkþ1 ; kþ1 1þ0:5sa2 V1;j kþ1 N1;j we get
1 1 kþ1=2 kþ1 kþ1 ¼ Gkþ1 ; bkþ1 ¼ bkþ1 : Gkþ1 akþ1 i;j i;j1 þ Fi;j1 i;j i;j1 ai;j1 i;j1 ai;j1 Find the first coefficients from (19) at j ¼ 1 and from the condition cz jH=2 ¼ 0 : 1 ~ kþ1 ; bkþ1 ¼ Fkþ1=2 =2: 1þV akþ1 i;1 1;1 i;1 ¼ 1 kð2Dkþ1 Þ i;1 The immediate task is to find the values of ckþ1 i;N2 , i ¼ 1,., N1 1 which is necessary to launch the sweep method. Write down the approximation of the boundary condition (4) ðz ¼ HÞ : k 0; . h
2
kþ1=2 kþ1=2 ð0:5sÞ ¼ 0:5 bkþ1=2 q2i qkþ1 2i q2i i 2 Dkþ1=2 cz tkþ1=2 ; ri ; H bkþ1 qkþ1 Dkþ1 cz ðtkþ1 ; ri ; HÞ : 2i
ð21Þ
In the boundary node with a precision of up to Oðh2z Þ we have
kþ1=2 kþ1 kþ1 kþ1 ckþ1 ¼ 0; k 0: 1;j1 G1;j c1;j þ c1;jþ1 þ F1;j
(16) th
Find approximation of concentration on the (k þ 1) layer in the form
time
(17)
The sweep coefficients are: j ¼ 2; 3; .; N2 ; akþ1 1;j
1 ~ kþ1 ; bkþ1 ¼ Fkþ1=2 =2: 1 þ V ¼ 1 kð2D Þ akþ1 kþ1 1;1 1;1 1;1 1;1
kþ1 kþ1 ckþ1 ¼ akþ1 i;j i;jþ1 ci;jþ1 þ bi;jþ1 : j ¼ 0; 1; .; N2 1; k 0:
~ kþ1 ckþ1 1W akþ1 1;j 1 1;j
kþ1 kþ1 kþ1 ckþ1 1;j ¼ a1;jþ1 c1;jþ1 þ b1;jþ1 ; j ¼ 0; 1; .; N2 1; k 0:
1243
1 1
kþ1=2 kþ1 kþ1 kþ1 Gkþ1 ¼ Gkþ1 ; bkþ1 : 1;j ¼ b1;j1 þ F1;j1 1;j1 a1;j1 1;j1 a1;j1
Find the first coefficients from (16) at j ¼ 1 and from the condition cz jH=2 ¼ 0 :
m m 2hz cz ðtm ; ri ; HÞzcm i;N2 2 4ci;N2 1 þ 3ci;N2 ; m ¼ k þ 1=2; k þ 1:
(22)
Concentration values on the (k þ 1/2)th layer are already known. For the (k þ 1)th layer, substituting values ckþ1 i;N2 2 and from expression (17) (at i ¼ 1) and (20) (at i ¼ 2,., N1 1), ckþ1 i;N2 1 we get the approximation
1 h kþ1 kþ1 3 þ akþ1 cz ðtkþ1 ; ri ; HÞz ckþ1 i;N2 ai;N2 1 4 i;N2 þ bi;N2 1 2hz
i kþ1 þ akþ1 i;N2 1 4 bi;N2 : Write down the approximation in a more compact form:
1244
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 6 ( 2 0 1 1 ) 1 2 3 9 e1 2 4 7
Fig. 2 e The impact of the parameter b0 .
Fig. 1 e Diffusion flux approximations of IðtÞ.
kþ1 kþ1 kþ1 kþ1 Ah3 þ akþ1 i;N2 ai;N2 1 4 ; Bhbi;N2 1 þ ai;N2 1 4 bi;N2 ; 2hz cz ðtkþ1 ; ri ; HÞzAckþ1 i;N2 þ B:
ð23Þ
Expressions (22), (23) substitute into (21), denote skþ1 i;N2 ¼ y: bkþ1 2 4 Dkþ1 A y þ G ¼ 0; y þ þ sgkþ1 2hz g2kþ1 ! bkþ1=2 kþ1=2 2 3Dkþ1=2 4 c ckþ1=2 Gh 2 þ 2hz sgkþ1=2 i;N2 gkþ1=2 i;N2
1 kþ1=2 kþ1=2 þ Dkþ1 B þ Dkþ1=2 ci;N2 2 4ci;N2 1 : 2hz
1 m qm 2i ¼ gm ci;N2 ; m ¼ k þ 1=2; k þ 1;
Taking the physical meaning of the process into account we choose the positive root of the square equation with respect to y. The roots have different signs, at least at relatively small s (denote s 1). The approximation order of the boundary condition is Oðh2z þ s2 Þ. It agrees with the approximation order of the scheme in the bulk. Now find all values of ckþ1 i;j at i ¼ 0 and i ¼ N1, j ¼ 0, 1,., N2. Using boundary condition on the axis of the cylinder kþ1 kþ1 kþ1 ðcr jþ0 ¼ 0Þ, we get ckþ1 0;j ¼ ð4c1;j c2;j Þ=3. Values of cN1 ;j > 0 are defined from the equation
h i n kþ1=2 kþ1=2 kþ1=2 kþ1=2 ci;j wkþ1 ¼ 4s1 wi;j þ a1kþ1=2 1 Wi;j a2kþ1=2 wi;j i;j ~ kþ1 ckþ1 akþ1 þ 4s1 1 : 1W þ akþ1 1 2 i;j i;j ~ s (s ¼ k þ 1/2; k þ 1) is the following: we The correction of W i;j s ~ ¼ ws =wmax and proceed computing using the assign W i;j i;j scheme of the method of alternating directions until stabili~ s zWs (usually twoethree iterations). zation W i;j i;j
3.4.
Classical diffusion problem
For simplified (rougher) calculations one can suppose that diffusion is the only limiting process (physicalechemical processes on the surface are ignored); then the Dirichlet boundary-value problem with zero boundary concentrations is considered:
v2 c 1 vc v2 c vc ; r˛ð0;LÞ; z˛ð0;HÞ; t˛ð0; t Þ; þ $ þ ¼ DðtÞ vt vr2 r vr vz2 cr ðt; þ0; zÞ ¼ 0; cðt; L;zÞ ¼ 0; cðt; r;0Þ ¼ cðt; r; HÞ ¼ 0; cð0;r; zÞ ¼ c: Partial sum of series
bkþ1 kþ1 2 2 3Dkþ1 kþ1 Dkþ1 kþ1 cN1 ;j þ cN1 2;j 4ckþ1 cN1 ;j þ þ N1 1;j 2 2hr 2hr sgkþ1 gkþ1 2 kþ1=2 c ¼ 0; s gkþ1=2 N1 ;j which approximates the condition (3) at r ¼ L. Using ckþ1 i;j compute
Table 1 e Parameter values common for all curves.a b0 ¼ 1013 D0 ¼ 106 g0 ¼ 103
Eb ¼ 12 104 ED ¼ 5 104 Eg ¼ 0
T0 ¼ 293 Tcrit ¼ 700 Te ¼ 773
a1 ¼ 0:1 a2 ¼ 0:2 a3 ¼ 103
L ¼ 4 103 H ¼ 2 102 c ¼ 1023
a ½b0 ; D0 ¼ m2 s1 , ½Eb ; ED ; Eg ¼ J mol1 , ½g0 ¼ m1 , ½T0 ; Tcrit ; Te ¼ K, ½a1 ; a2 ; a3 ¼ s1 , ½L; H ¼ m, ½c ¼ m3 .
Fig. 3 e The impact of the parameter Eb .
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 6 ( 2 0 1 1 ) 1 2 3 9 e1 2 4 7
Fig. 4 e The impact of the parameter D0 .
8 9 Zt < = m0
np cðt; r; zÞ ¼ Anm exp lnm DðsÞds $J0 m r $sin z ; : ; L H n¼1 m¼1 0
np2 m0 2 4 1 þ ð1Þnþ1 c ; lnm h þ m Anm h L H pnm0m J1 m0m N X N X
is used as approximation of the solution. This sum is a generalized solution since boundary conditions are inconsistent. Here, DðtÞhDðTðtÞÞ, J0, J1 are zero-order and first-order Bessel functions of the first kind, respectively, m0m are sequential zeros of the function J0(m). H flux through the surface is: Z Dcn ds ¼ c
IðtÞ ¼
i h nþ1 N X N 8HDðtÞ 1 þ ð1Þ X
np n¼1 m¼1 " 2 # 1 þ ð1Þnþ1 pL $expf.g: þ ð1Þnþ1 2n n Hm0m
1245
Fig. 6 e The impact of the parameter g0 .
time is not so important: interesting effects take place in the neighborhood of the TDS-peak “top”. How many terms of a series should we include in the calculations for appropriate curve approximation (5%)? The series converges slowly. For instance, at D0 ¼ 2 103 m2 s1, ED ¼ 6 104 J mol1 (Dz2 107 , T ¼ 773 K ), c ¼ 1023 m3 one should take 144 items of a series (n ¼ 1, 3, 5, 7, 9, 11, 13, 15; 1 m 18). The series equals zero when n is even. Fig. 1 shows the approximation of I(t) with partial sums (n, m) (N, M ) ¼ (15, 18), (5 10), (3, 3), (1, 1) with respect to the decrease of the maximum. It is insufficient to include several items, the peak is unique (in the concrete situation it is necessary to take several comparable diffusion channels with different D( j ) into consideration).
S
4. Algorithm and results of numerical modelling
One often takes only the first item (n ¼ m ¼ 1) into account, believing that the rest of the items are quickly damped. Yet, in the problem under consideration the asymptotic form versus
Describe the computing algorithm step-wise. We define the values of L, H, D0, ED, b0, Eb, g0, Eg, c, wmax, g, T0, Te, Tcrit, a1, a2,
Fig. 5 e The impact of the parameter ED .
Fig. 7 e The impact of the parameter Eg .
1246
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 6 ( 2 0 1 1 ) 1 2 3 9 e1 2 4 7
Fig. 9 e The impact of the defects of two types. Fig. 8 e Comparison of JðtÞ and IðtÞ.
a3, s, r, ~c. Transition from the kth layer to the (k þ 1)th layer is realized in two steps. Step I: the computing algorithm on the (k þ 1/2)th layer ðk 0Þ is the following. ~ from equation (9) at 1. Compute the values of W i;j s ¼ k þ 1=2. 2. According to expressions (10), (11) compute the sets of kþ1=2 kþ1=2 , bi;j . coefficients ai;j 3. Determine the values of concentration in boundary nodes, kþ1=2 solving the square equation with respect to y ¼ sN1 ;j > 0. 4. Using formula (11) find the approximate values of concentration in all internal nodes. kþ1=2 in 5. Using boundary conditions determine the values of ci;j boundary nodes at j ¼ 0 and j ¼ N2 , i ¼ 0; 1; .; N1 . 6. From the second equation of the difference scheme in the bulk compute values of concentration in defects of the reversible capture type:
6. From the second equation of the difference scheme in the bulk obtain: h i N kþ1=2 kþ1=2 kþ1=2 kþ1 wi;j ¼ kþ1 ; N ¼ 4s1 wi;j þ a1kþ1=2 1 Wi;j ci;j a2 þ 4s1 kþ1=2 ~ kþ1 ckþ1 : 1W a2kþ1=2 wi;j þ akþ1 i;j 1 i;j
kþ1=2
kþ1=2 ¼ wi;j
h i kþ1=2 k ~ kþ1=2 ckþ1=2 1 W 4s1 wki;j þak1 1Wi;j cki;j ak2 wki;j þa1 i;j i;j
7. Correction
a2kþ1=2 þ4s1 of
~ kþ1=2 W i;j
is
the
following:
one
:
assign
~ kþ1=2 ¼ wkþ1=2 =wmax and proceed as in items 2e6 until W i;j i;j ~ kþ1=2 zWkþ1=2 (twoethree iterations). stabilization W i;j i;j Step II: the computing algorithm on the (k þ 1)th time layer ðk 0Þ is the following. ~ from equation (9) at s ¼ k þ 1. 1. Compute the values of W i;j 2. According to expressions (16), (17), (19), (20) compute coefkþ1
kþ1 kþ1 , bkþ1 ficients akþ1 1;j , ai;j , bi;j . 1;j
3. Determine the values of concentration in boundary nodes, solving the square equation with respect to y ¼ ckþ1 i;N2 > 0. 4. Using formulas (17), (20) find the approximate values of concentration in internal nodes. at i ¼ 0 and 5. From boundary conditions determine ckþ1 i;j i ¼ N1 , j ¼ 0; 1; .; N2 .
7. Correction
of
~ kþ1 W i;j
is
the
following:
one
assign
~ kþ1 ¼ wkþ1 =wmax and proceed as in items 2e6 until stabiW i;j i;j kþ1 ~ lization W i;j zWi;j . kþ1
4.1.
Results of numerical modelling
Desorption flux curves JðtÞ are represented in the figures. Varied coefficient values correspond to the curves in the following way: the first value is for the curve with the highest maximum, the last one is for that with the lowest maximum. The small circle shows the time instant when the initial surface hydrogen is exhausted. The algorithm developed for numerical modelling of hydrogen thermodesorption is applicable in a wide range of parameters. In our figures however, we limit ourselves to concrete values of the parameters presented below. We define g ¼ 1023 m3 , wmax ¼ 1023 m3 . The default parameters, which are common for all curves, are taken from the Table 1. Fig. 2 illustrates the influence of the coefficient b0 on the desorption flux curve. Changes are: b0 ¼ 1012 ; 1013 ; 1014 m2 s1 . The influence of the activation energy Eb can be seen in Fig. 3, changes in order of decrease of the maximum are: Eb 104 ¼ 9; 10; 12 J mol1 . Changes for Fig. 4 are: D0 ¼ 5 105 ; 105 ; 106 m2 s1 . Fig. 5 shows the influence of the parameter ED : ED 104 ¼ 4; 5; 6. The variations of g0 ¼ 50; 100; 200 m1 are seen in Fig. 6. Fig. 7 shows the dependence of the desorption curve on the value of Eg 103 ¼ 1; 0; 1. Fig. 8 shows two flux curves: the left-hand curve is the desorption flux for the accepted boundary-value problem with dynamic boundary conditions, the right-hand curve is the diffusion flux for the degenerate Dirichlet problem. Parameter values are: D0 ¼ 2 103 , ED ¼ 6 104 , b0 ¼ 105 , Eb ¼ 12 104 , g0 ¼ 103 , Eg ¼ 0, a10 ¼ 0, a3 h0, g ¼ 3 1023 m3 ,
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 6 ( 2 0 1 1 ) 1 2 3 9 e1 2 4 7
wmax ¼ 1024 m3 . The flux curve taking the defects of two types into account is illustrated in Fig. 9 (reversible capture and decomposition).
Acknowledgements The work was supported by the Russian Foundation for Basic Research (grant 09-01-00439) and by the Russian Science Support Foundation.
[4]
[5]
[6]
references
[1] Kolachev BA. Hydrogen embrittlement of metals. Moscow: Metallurgy; 1985 [in Russian]. [2] Kunin LL, Golovin AI, Surovoi UI, Hohrin VM. Problems of metals degassing. Moscow: Nauka; 1972 [in Russian]. [3] Polyanskiy AM, Polyanskiy VA, Popov-Diumin DB, Kozlov EA. New measuring complex for absolute determination of
[7]
[8]
1247
hydrogen content in hydrogen energetics materials. Int J Altern Energy Ecology 2006;6:29e31. Gabis IE, Kompaniets TN, Kurdyumov AA. In: Zakharov AP, editor. Interactions of hydrogen with metals. Moscow: Nauka; 1987. p. 177e206 [in Russian]. Polyanskiy AM, Polyanskiy VA, Yakovlev Yu A. Methodic of investigation of the binding energies of hydrogen in solid state on the basis of hydrogen analyzer AV-1. Proc. III Int. Conf. “Interaction of hydrogen isotopes with structural materials”. Sarov, RFNCeVNIIEF; 2008. p. 346e53. Zaika YuV, Rodchenkova NI. Modelling of diffusion TDSspectrum peak of dehydriding with size reduction and heat absorption effects. In: Baranowski B, et al., editors. Carbon nanomaterials in clean energy hydrogen systems, NATO science for peace and security series e C: environmental security. Springer; 2008. p. 863e78. Zaika YuV, Rodchenkova NI. Boundary-value problem with moving bounds and dynamic boundary conditions: diffusion peak of TDS-spectrum of dehydriding. Appl Math Model 2009; 33:3776e91. Samarskij AA. Introduction to difference schemes theory. Moscow: Nauka; 1971 [in Russian].