# The Inverse Problem of Identification of Hydrogen Permeability Model.

1. IntroductionStudies on the interaction of hydrogen isotopes with structural materials are mainly necessitated by problems in the energy industry, metal protection from hydrogen corrosion, and the design of chemical reactors [1-10]. The limiting factors are diffusion processes as well as physical and chemical phenomena at the surface. The transfer parameters depend on the technological features of material batch production. It is therefore unreasonable to target at "tabular data". Instead, effective algorithms for solving the inverse problems of parametric identification of adequate mathematical models by experimental curves (data) are necessary. In this study, we consider the permeability model taking into account the main factors and the self-descriptiveness of the experiment. We shall focus on the methods of permeability and thermal desorption taking into account only the main limiting factors for the applied membrane filtering problem and the informative capabilities of the considered experimental methods. The mathematical research is based on the articles [11-13], which provide descriptions of the experimental techniques and experimental material on promising alloys for hydrogen separation and purification.

One had to estimate the parameters of diffusion and sorption to numerically model the different scenarios and experimental conditions of the material usage (including extreme ones) and identify the limiting factors. Some particular problems of the hydrogen materials science related to the topic of this study were presented and investigated in [14-23]. This work is a continuation of [24-27], where the results of modelling hydrogen thermal desorption under various limiting factors are presented. This article deals with the inverse problem of parametric identification based on the suggested "cascade" experiment.

Experimental practices usually employ various modifications of the penetration method and TDS. The results of measurements depend both on the unit design features and on the procedure of preparing samples for hydrogen permeability testing. A successive use of various methods often causes, for example, impurities to appear on the sample surface, which significantly affects the reproducibility of the results. These data are the input for the inverse problems of parametric identification, which are sensitive to the level of different errors. It is therefore advisable to aggregate experiments to improve the accuracy and informative value of the measurements. We suggest the following set-up of the "cascade" experiment.

A membrane heated to a fixed temperature served as the partition in the vacuum chamber. Degassing was performed in advance. A sufficiently high pressure of hydrogen gas was built up in steps at the inlet side. The penetrating flux was determined by mass-spectrometry in the vacuum maintained at the outlet side. This is the penetration method. Its advantage is a reliable determination of the diffusion coefficient by the Daines-Berrer method (based on the so-called lag time). It allows distinguishing between the bulk and the surface processes in the model, keeping in mind that surface parameters are significantly more difficult to estimate. When the steady state level of the penetrating flux is registered, we increase the inlet pressure and wait until a new steady state value is established. Using (at least) three pressure jumps at the inlet side we record the steady state flux values at the outlet side, thus determining "the degree of rectilinearity" of the isotherm. Then the pumping for vacuum building is stopped and the experiment proceeds as the "communicating vessels" method. When pressure values become nearly equal (the sample is almost uniformly saturated with hydrogen) it is possible to turn off the heating, create the vacuum at both sides of the membrane, and begin slowly reheating the sample (TDS-experiment). In addition, there is no depressurization of the diffusion cell and the sample surface remains uncontaminated by additional impurities. We will clarify the details of the aggregated experiment stages as we describe the method of solving the inverse problem of parametric identification. An important consideration is the uniqueness of the parameter estimates of the investigated model. Mathematicians are often reproached for "fascination with uniqueness theorems". But after all, in justifying the choice of, for example, structural materials for the ITER project, the results obtained on thin laboratory samples are extrapolated to "walls". Uniqueness allows for a correct recomputation.

Papers [18-20, 27] were dedicated to interpretation of TDS peaks. Analysis of the causes of various TDS surges is crucial for the selection of reactor structural materials contacting with hydrogen isotopes. A sufficiently detailed review is presented in [19, 20]. Where modeling does not take the dynamics of surface processes into account, TDS peaks are inevitably interpreted as a result of capture of diffusing atomic hydrogen by the structural defects (traps) of the material with different binding energies and (or) as a manifestation of multichannel diffusion [28]. The above listing of causes is not exhaustive. In this article, modelling shows that the appearance of a desorption peak can be caused by some combinations of the rates of surface and diffusion processes. This complicates the problem of TDS spectrum interpretation even more.

The main result of the paper is the method of parametric identification from experimental data. The difficulties of inverse problems solving in mathematical physics are well known. There is extensive mathematical literature and a number of specialized journals (inverse problems, ill-posed problems, etc.). In experimental practice, the inverse problem of multiparameter estimation is reduced to the one-factor-at-a-time method for DLR and SLR. In real life, however, a material is used in the presence of a dynamic "surface-bulk" interplay. Thoroughly elaborated techniques are available for estimation of the diffusion coefficient. The determination of desorption and dissolution parameters is far more complex (unless the temperature is artificially lowered to "turn off" processes in the bulk). The paper presents an algorithm allowing the estimation of desorption and dissolution coefficients where diffusion and surface processes interact intensively.

2. Hydrogen Permeability Model

2.1. Distributed Transfer Model. Let us briefly describe the experiment. A sample of a structural material preheated to a fixed temperature acts as a vacuum chamber barrier. The sample degassing is performed in advance. At the initial time moment, pressure is built up at the inlet side by puffing of a portion of molecular hydrogen. The declining pressure in the input chamber and increasing pressure in the output chamber are measured.

Consider hydrogen transfer through the sample (l is the plate thickness and S is its area). The temperature T is constant throughout the experiment. The concentration of dissolved diffusing hydrogen (in monatomic state) is sufficiently low and the diffusion flux can be considered proportional to the concentration gradient. The membrane is thin and the material has a sufficiently high hydrogen permeability coefficient, so we restrict ourselves to a standard diffusion equation:

[[partial derivative].sub.t]c(t,x) = D(T)[[partial derivative].sup.2.sub.x]c(t,x), t > 0, c(0,x) = 0, x [member of] [0,l], (1)

where c(t, x) is the concentration of diffusing hydrogen. The diffusion coefficient D depends on the sample temperature T in an Arrhenius way with preexponential factor [D.sub.0] and activation energy [E.sub.D]: D = [D.sub.0] exp[-[E.sub.D]/[RT(t)]}.

Initial data are determined by the fact that the sample had been preliminarily degassed: c(0, x) = 0, z(0, x) = 0, x [member of] [0,l].

Nonlinear boundary conditions are derived from the material flux balance:

[mathematical expression not reproducible]. (2)

Here, [Q.sub.in](i), [Q.sub.out](t) are the amounts of hydrogen atoms in the input chamber of volume [V.sub.in] and output chamber of volume [V.sub.out], pc.sub.0](t) [equivalent to] c(t, 0), [c.sub.l](t) [equivalent to] (t,l). The identity sign is frequently used here in the sense of equality by definition. Within the considered operating temperature range the gaseous hydrogen is in molecular form, but for consistency (considering that atomic hydrogen diffuses through the metal membrane) we use atoms as the unit. According to the kinetic gas theory, the incident particle flux density [J.sub.p] is related to the pressure p by the Hertz Knudsen formula: [J.sub.p] = p/ [square root of 2[pi]mkT] (k is the Boltzmann constant, m is the mass of a hydrogen molecule). In the context of the experiment it is convenient to choose the following measurement units [l] = cm, [p] = Torr. Then we numerically obtain the dependence [mathematical expression not reproducible]. The processes of physical adsorption, chemisorption, dissociation of molecules into atoms, and dissolution take place on the surface. Only a small part of "incident" H atoms will, however, be absorbed into the membrane volume. This is taken into account by the factor 2s. One can write s (as a parameter of the model) instead of 2s, but it is more convenient to interpret the dimensionless probability factor s as the fraction of absorbed H atoms within the 2s notation. Thus, 2s[mu]p is the resulting flux of atoms through the surface into the bulk without differentiation into more elementary stages. We will omit the word "density" assuming that the surface has unit area.

Hereinafter, [J.sub.0,l] = b[c.sup.2.sub.0,l] are the densities of the desorption flux from the sample (deviation from the square desorption law is significant only at extreme temperatures) and b is the desorption coefficient. We also assume that s and b depend on the temperature in an Arrhenius way. Formally, the activation energy [E.sub.s] in the exponent can as well be negative, being a linear combination of the activation energies and heats of the surface processes on the way "from gas to the solution". If a constant saturation pressure of molecular hydrogen [p.sub.s] = const is maintained at a constant temperature T = const on both sides of the membrane, the equilibrium concentration [bar.c] of the dissolved atomic diffusionally mobile hydrogen is finally established. By equating the derivatives in model (2) to zero, we get Sieverts' law [bar.c] [varies] [square root of [p.sub.s]] : [bar.c] = [GAMMA] [square root [p.sub.s]], [GAMMA] [equivalent to] [square root of 2s[mu]/b].

Let us clarify the experimental conditions. The volumes [V.sub.in,out] comprise several liters, the thickness of the membrane is l < mm, the area S is about 1 [cm.sup.2], and the inflow pressure [p.sub.0](0) is within several hundreds Torr.

It now remains to find the magnitudes of [Q.sub.in], [Q.sub.out]. Within the time of transfer through the membrane the gas is in the thermodynamical quasi-equilibrium, wherefore we use the formula N = pV/(kT). Here, N is the number of gas particles occupying the volume V at the temperature T and the pressure p (in the SI system [p] = Pa, [V] = [m.sup.3], [k] = J/K). Taking into account the relationships Torr = 133.322 Pa and Pa = J/[m.sup.3] (formally), we get the following expressions for the corresponding pressures and volumes in the boundary conditions (2): Q = 2N = [alpha]pV/T, [alpha] [approximately equal to] 1.931 * [10.sup.19]. Here, p, V, T are the numerical values in the selected units (Torr, [cm.sup.3], K).

Within the experimental unit the membrane is situated in a tube (which is heated to a predetermined temperature) between the inlet and the outlet chambers. The tube diameter is large enough to consider the equality of pressures as the criterion of thermodynamical quasi-equilibrium between the gas in the tube and in the chambers. The membrane temperature should be taken for the formula for the kinetic constant [mu](T). The gas inside the volumes [V.sub.in,out] (whose massive walls are at room temperature) may get heated up. During the preliminary experiment it is recommended to fill the chambers with a practically impermeable metal membrane between them with "ambient" gas. Then we heat the tube and record the pressure rises. Within the framework of the ideal gas approximation (equation of state) this procedure enables estimation of the increments of the gas temperature inside the chambers. The corresponding gas temperatures are the ones to be used in the formula given for Q (and the subsequent ones, only excluding the value [mu]). The need of such a refinement arises from the characteristics of this particular experimental unit. Such an adjustment of the values of T should not cause difficulties in further calculations. Besides, this procedure has relatively little effect on the final calculation of the model pressures taking into account the measurement errors and relatively large volumes V.

2.2. Fast Hydrogen Permeability Model. It is clear from physical considerations that a quasi-stationary state is quickly established when the membrane is thin and the material has a sufficiently high hydrogen permeability coefficient: the diffusant concentration distribution is practically linear with respect to the thickness. In this sense, the results of numerical modelling based on the "general" model (the presented boundary-value problem) confirm its adequacy. Since near-to-surface concentrations cannot be measured, the Richardson approximation is usually used in practice to analyze the penetrating flux:

J(t) = -D[[partial derivative].sub.x]c = D[l.sup.-1][[c.sub.0](t) - [c.sub.l](t)] [approximately equal to] [J.sub.R](t) = D[GAMMA][l.sup.-1][[square root of [p.sub.0](t)] - [square root of [p.sub.l](t)]]. (3)

Let us formulate the problem of modelling the concentrations [c.sub.0,l] using the pressures [p.sub.0,l] (the problem is also of interest per se) without the quasi-equilibrium simplification c(t) = [GAMMA] [square root of p(t)]. The quasi-stationary state is achieved within a time [t.sub.0], which is short compared to the total experiment time ([[partial derivative].sub.x]c = -[[c.sub.0](t) - [c.sub.l](t)]/l). So, the original model (1)-(2) can be simplified (taking into account the formula Q = [alpha]pV/T):

[mathematical expression not reproducible], (4)

2s[mu][p.sub.0,l](t) - b[c.sup.2.sub.0,l](t) = [+ or -] D[l.sup.-l][[c.sub.0](t) - [c.sub.l](t)], t [greater than or equal to] [t.sub.0] > 0. (5)

Since by virtue of the "inlet-outlet" balance the equalities hold true,

[[??].sub.l](t) = -[V.sub.in][V.sup.-1.sub.out][[??].sub.0](t) [??] [p.sub.l](t) = [p.sub.l]([t.sub.0]) + [V.sub.in][V.sup.-1.sub.out][[p.sub.0]([t.sub.0]) - [p.sub.0](t)], (6)

it is sufficient to express the concentrations [c.sub.0,l](t) = [c.sub.0,l]([p.sub.0](t)) from the boundary conditions (5) and substitute them into the first equation of (4) (the sign is chosen depending on whether the index is 0 or l). The dimensionless normalized variables are convenient for numerical simulation:

[X.sub.0,l](t) = 1 + 2l[c.sub.0,l](t)b[D.sup.-1], [a.sub.0,l](t) = 4[l.sup.2][[GAMMA].sup.2][p.sub.0,l](t)[b.sup.2][D.sup.-2] - 1. (7)

In addition, the system of equations (5) is compactly written in the symmetric form [a.sub.0] + 2[X.sub.l] = [X.sup.2.sub.0], [a.sub.l] + 2[X.sub.0] = [X.sup.2.sub.l]. For the variable X [equivalent to] [X.sub.l] we obtain the incomplete quartic equation [[[X.sup.2] - [a.sub.l]].sup.2] = 4[2X + [a.sub.0]], which can be solved in radicals (for physical considerations we are interested in the positive root). However, the explicit expression is somewhat cumbersome and we will anyway have to numerically integrate the first equation of (4) in the form [[??].sub.0] = f([p.sub.0]). Therefore, we shall aim to derive differential equations for [X.sub.0,l], since information about the dynamics of the boundary concentrations [c.sub.0,l] is also important.

Differentiate (5) with respect to time and substitute the pressure derivatives from (4). For the variables [X.sub.0,l] we get the system

[[??].sub.0](t) = -2s[M.sub.0][[X.sub.0] - [X.sub.l]] x [X.sub.l] - [V.sub.in][V.sup.-1.sub.out]/[X.sub.0][X.sub.l] - 1, [M.sub.0] [equivalent to] [mu]sT/[alpha][V.sub.in], (8)

[[??].sub.l](t) = -2s[M.sub.l][[X.sub.0] - [X.sub.l]] x [X.sub.0] - [V.sup.-1.sub.in][V.sub.out]/[X.sub.0][X.sub.l] - 1, [M.sub.l] [equivalent to] [mu]sT/[alpha][V.sub.out]. (9)

The change of variables in (7) determines the concentrations [c.sub.0,l](t) from which the model pressures [p.sub.0,l](t), t [greater than or equal to] [t.sub.0] are calculated using (5).

Let us formulate step-by-step the numerical algorithm of modelling the pressures [p.sub.0,l](t) for the current values of D, b, s coefficients (the authors used the Scilab freeware). We target at the "normal" experimental conditions [11,12,29-31], including the values of p, T, l, V, S.

(1) We fix t = [t.sub.0]: omit fast transient processes (the duration of the transient processes is about tens of seconds on the hours-long experimental time scale). For the variable X [equivalent to] [X.sub.l] we choose the root of the biquadratic polynomial [[[X.sup.2] - [a.sub.l]([t.sub.0])].sup.2] - 4[2X + [a.sub.0]([t.sub.0])]. From physical considerations it follows that [c.sub.l]([t.sub.0]) > [[bar.c].sub.l]([t.sub.0]) and thus X > 1 + [square root of 1 + [a.sub.l]] = 1 + 2l[GAMMA] [square root of [p.sub.l]([t.sub.0])b[D.sup.-1].

(2) The system of equations [a.sub.0] + 2[X.sub.l] = [X.sup.2.sub.0],[a.sub.l] + 2[X.sub.0] = [X.sup.2.sub.l](t = [t.sub.0]) yields the missing value of [X.sub.0]([t.sub.0]). Formally, one equation is enough, but we take into account averaging procedures including determination of the values of [p.sub.0,l]([t.sub.0]).

(3) We numerically integrate the ODE system (8), (9) with the obtained initial data. The change of variables in (7) defines the concentrations [c.sub.0,l](t), which are used to calculate the model pressures [p.sub.0,l](t), t [greater than or equal to] [t.sub.0] from (5).

Computational experiments show that the model curves almost coincide (at t [greater than or equal to] [t.sub.0]) with those generated by the originally proposed model, i.e., the nonlinear distributed initial boundary value problem.

Observe the difference from the quasi-equilibrium model (the Richardson approximation), where the only parameter for approximation of the experimental pressure is the complex [PHI] = D[GAMMA]. All the variable parameters of the original model that influence the permeability, namely, D, b, s, are important when running the above algorithm. Thus, the fast hydrogen permeability model does not lose in informativeness concerning the considered transfer parameters.

3. Modelling of Hydrogen Permeability

3.1. Numerical Modelling of the Penetration Experiment. The proposed model is adapted to the experimental conditions and the data range for alloys based on V group metals with high hydrogen permeability, in particular, data for vanadium alloys which are presented in [11-13,29-33]. We fix T = 673 K, [mathematical expression not reproducible], D = 2 x [10.sup-5] [cm.sup.2]/s, [GAMMA] = 2 x [10.sup.20] [1.sub.H]/([cm.sup.3] [square root of Torr]), [PHI] = D[GAMMA] = 4 x [10.sup.15] [1.sub.H]/(cms [square root of Torr]). We set the value s = 1.2 x [10.sup.-4] and calculate the corresponding desorption coefficient b = 2[mu]s/[[GAMMA].sup.2] = 5.7 x [10.sup.-24] [cm.sup.4]/s. The input pressures [[bar.p].sub.1,2,3] = {30, 50,70} Torr were built up in steps at the inlet and maintained to achieve steady state fluxes at the outlet. Degassing of the membrane was done in advance and continuous vacuum pumping was performed at the outlet side. The gas temperature inside the inlet and outlet chambers (whose volumes are sufficiently high) is assumed to be equal to 300 K. This slight difference from the room temperature is due to the heating of the diffusion cell with the sample inside (specified by the characteristics of the equipment). The experimental conditions are such that the concentration at the membrane outlet side is near zero and at the inlet side a stationary concentration is quickly established (but it is lower than the equilibrium one): [??] < [bar.c]. Within the model we determine [[bar.c].sub.i] and [[??].sub.i] by the formulas

[mathematical expression not reproducible]. (10)

For the given values of the parameters we have [[??].sub.1] = 1.06 x [10.sup.21] < [[bar.c].sub.1] = 1.09 x [10.sup.21], [[??].sub.2] = 1.38 x [10.sup.21] < [[bar.c].sub.2] = 1.41 x [10.sup.21], [[??].sub.3] = 1.64 x [10.sup.21] < [[bar.c].sub.3] = 1.67 x [10.sup.21].

Next we express the solutions of the standard boundary value problems with Dirichlet boundary conditions corresponding to the jumps of the inlet pressure.

Stage I. The boundary value problem of the penetration method is

[mathematical expression not reproducible]. (11)

The penetrating flux is [J.sub.1] = -D[c.sub.x|l] = D[[??].sub.-1][l.sup.-1][1 + 2 [[summation].sup.[infinity].sub.n=1][(-1).sup.n] exp{-[n.sup.2][[pi].sup.2]Dt/[l.sup.2]}].

From the computational point of view it is convenient to introduce the dimensionless time t' = Dt/[l.sup.2] which is oriented at the characteristic diffusion time [l.sup.2]/D. As small t [right arrow] 0 a singularity appears if we directly use the partial sum of the expression

f(t') = 1 + 2 [[infinity].summation over (n=1)] [(-1).sup.n] exp {-[n.sup.2][[pi].sup.2]t'}. (12)

Let us provide another expression for f using the properties of the Jacobi theta function. More precisely, we are interested in the function

[[theta].sub.3](t,x) = 1 + 2 [[infinity].summation over (n=1)] exp {-[n.sup.2][[pi].sup.2]t} cos (2n[pi]x), t > 0. (13)

We have an alternative presentation for x = 0 [34]:

[mathematical expression not reproducible]. (14)

The series on the left is rapidly converging for large t. But the series on the right is rapidly converging for small t. If we define [theta](t) = [summation] exp{-[pi][n.sup.2]t} (n [member of] Z, t > 0), then we get [theta](1/t) = [square root of t][theta](t), or [square root of t] [summation] exp{-[pi][n.sup.2]t} = [summation] exp{- [pi][n.sup.2][t.sup.-1]} (n [member of] Z), which is known as the functional equation for the theta function [35]. After some auxiliary transformations for 0 < t' << 1 we get the appropriate representation

f(t') = 2/[square root of [pi]t'][summation over (m=2n-1)] exp{-[m.sup.2]/4t'} (n [member of] N). (15)

The function f (t') has an S-shaped form of the saturation curve (see the inset in Figure 1). Stage I ends with c([t.sub.*] x) = [[??].sub.1](l - x)[l.sup.-1].

Stage II. [t.sub.*] [right arrow] [t.sub.0] = 0 is the new t time zero:

[mathematical expression not reproducible]. (16)

Stage II ends with c([t.sub.*],x) = [[??].sub.2](l - x)[l.sup.-1]. The penetrating flux is

[mathematical expression not reproducible]. (17)

Stage III. The formulas are similar to the cyclic interchange [[??].sub.1] [right arrow] [[??].sub.2], [[??].sub.2] [right arrow] [[??].sub.3].

The result of connecting the stages into a single "experimental" curve of the penetrating flux J(t) (conventionally, t = ([t.sub.1],[t.sup.*.sub.1] + [t.sub.2], [t.sup.*.sub.2] + [t.sub.3]), J = ([J.sub.1],[J.sub.2], [J.sub.3])) is shown in Figure 1.

3.2. Converting Pressure into Flux. During the real experiment, the gas pressure inside the outlet volume V = [V.sub.out] is measured, not the penetrating flux. Therefore, in this subsection we will provide the corresponding recalculation formula. Denote the pumping rate of the vacuum system by v ([v] = [m.sup.3]/s). We take as the framework the ideal gas state equation: pV = NkT. Here, [p] = Pa, [V] = [m.sup.3], N is the number of particles ([H.sub.2] molecules), and k is the Boltzmann constant. Differentiating on t we get [??]V = NkT. Let us calculate the particle balance over the time [DELTA]t:

[DELTA]N = 0.5[J.sub.l]S[DELTA]t - v[DELTA]tN[V.sup.-1], [J.sub.l] = -SD[[partial derivative].sub.x]c (t, l), (18)

where [mathematical expression not reproducible] (the concentration), [[J.sub.l]] = [1.sub.H]/([m.sup.2] s). The factor 0.5 is due to the fact that the diffusion flux is atomic, and the particle in the volume V is the [H.sub.2] molecule. Let us divide the equation of the material balance by [DELTA]t and direct [DELTA]t [right arrow] 0. Finally, we get the differential equation:

[mathematical expression not reproducible], (19)

Denoting [[theta].sub.1] = kTS/(2V), [[theta].sub.0] = V/v, we get [??] = [[theta].sub.1][J.sub.l] - p/[[theta].sub.0] ([p.sub.0] = 0),wherefrom

[mathematical expression not reproducible]. (20)

The measured function p(t) is noised, so we first apply a smoothing procedure and only then calculate the derivative [??](t). Note that when the pumping system is powerful enough and permeability is relatively slow the first summand acts as a minor correction to the approximation [J.sub.l] [approximately equal to] 2vp(t)/(kTS). Asymptotically, a steady state value is established: p [bar.J] = 2v[bar.p]/(kTS).

Only the operation of integration is needed to calculate the lag time (see below):

[[integral].sup.t.sub.0] [J.sub.l]([tau])d[tau] = 2V/kTS p(t) + 2v/kTS [[integral].sup.t.sub.0] p([tau]) d[tau], (21)

so a preliminary approximation of the derivative [??](t) is unnecessary. It suffices to use the composite Simpson formula. One should remember that we use the following measurement units within this subsection [kT] = J = N x m, [L] = m, [p] = Pa = N/[m.sup.2] (Pa [approximately equal to] 7.5 x [10.sup.-3] Torr), [J] = [1.sub.H]/([m.sup.2] s). In the following, we return to the measurement units accepted in this article.

3.3. Boundary Value Problem of Hydrogen Permeability. When the steady state permeability value is established during the penetration experiment, continuous pumping at the outlet and maintenance of constant pressure at the inlet are stopped. The aggregated experiment moves to the stage of "communicating vessels": inlet pressure declines and outlet pressure grows ([p.sub.0,l](t) are measured). Considering the new time zero we also have t [greater than or equal to] 0. We are so far talking about the direct problem of modelling hydrogen pressures inside the volumes [V.sub.in,out]. We specify the values: S = 0.5 [cm.sup.2], [V.sub.in] = 1500 [cm.sup.3], [V.sub.out] = 2200 [cm.sup.3], [p.sub.0](0) = [[bar.p].sub.0] = [[bar.p].sub.3], [bar.c] = [GAMMA] [square root of [[bar.p].sub.0]]. We target here the experimental conditions and the data for the [V.sub.85][Ni.sub.15] alloy [11, 33].

We numerically solve the initial boundary problem of hydrogen permeability:

[mathematical expression not reproducible]. (22)

The membrane temperature is taken in the dependences of D, b, s, [mu] on T, and the temperature of the gas inside the chambers is taken in the expressions for [Q.sub.in,out] (take into account the correction to room temperature due to the heating of the diffusion cell). The model curves of molecular hydrogen pressures are presented in Figure 2. If we use the ODE system (8), (9) instead of the "full" model, standard software packages will suffice (we substitute the values of the hydrogen temperature inside the volumes [V.sub.in,out] into the expressions for [M.sub.0,l]). To this end one should skip the initial time [t.sub.0] within several minutes until a quasi-stationary (not quasi-equilibrium) mode is established. Then the above algorithm is applied to the fast hydrogen permeability model. Figure 2 shows that there is no significant additional error during the numerical simulation (the curves visually coincide). Model curves were numerically generated to test the following algorithm for solving the inverse problem of parametric identification. The parameters generating these curves were then "forgotten".

4. General Identification Technique

4.1. Determination of the Lag Time. For completeness, we briefly describe the method of estimating the diffusion coefficient proposed by Daines-Berrer. The curve of the flux [J.sub.1](t) asymptotically moves to the stationary value [mathematical expression not reproducible].

The intersection of the asymptote with the t axis gives the so-called lag time [[tau].sub.0] = [l.sup.2]/(6D), which allows estimating the diffusion coefficient. Analytically,

[mathematical expression not reproducible]. (23)

Note that the value under the integral sign is a relative magnitude which does not require absolute values ofthe penetrating flux in any measurement units ([[bar.J].sub.1] = sup [J.sub.1](t)). In addition, the value [[tau].sub.0] does not depend on [[??].sub.1]. It is usually assumed that the locally equilibrium concentration [[bar.c].sub.1] = [GAMMA] [square root of [[bar.p].sub.1]] is quickly established at the inlet, so one can additionally estimate the solubility [GAMMA] = [square root of 2s[mu]/b[ and the permeability [PHI] = D[GAMMA] using the corresponding value [[bar.J].sub.1] = D[[bar.c].sub.1]/l. This assumption is not used in this article. We weaken it to increase the accuracy of further estimations. We assume that according to the experimental conditions the stationary inlet concentration [c.sub.0](t) [approximately equal to] [[??].sub.1] < [[bar.c].sub.1] (t [greater than or equal to] [epsilon], [epsilon] << 1) is quickly established given that [c.sub.l](t) [approximately equal to] 0. The value of [[??].sub.1] as such is yet to be clarified. Thus, only the estimate of the diffusion coefficient D is considered reliable at this stage.

With a new zero time reference, integrating the expression [J.sub.2](t) we get

[[integral].sup.t.sub.0] [[J.sub.2]([tau]) - [[bar.J].sub.1]] d[tau] [approximately equal to] [[bar.J].sub.2] - [[bar.J].sub.1]] x [t - [l.sup.2][(6D).sup.-1]] ([t.sub.*] [right arrow] [t.sub.0] = 0), (24)

where [[bar.J].sub.i] = D[[??].sub.i]/l, t [greater than or equal to] [t.sub.*], = [t.sup.*.sub.2] Formally, we obtain the same expressions for the lag time and the estimate of D if we change both the zero time and the flux baseline ([[bar.J].sub.1] value overstatement). There is no additional information here (about the target values of the surface parameters b and s), but the triple penetration experiment allows refining D. For the model numerical experiment we have [[tau].sub.01] [approximately equal to] [[tau].sub.02] [approximately equal to] [[tau].sub.03] [approximately equal to] 20.82 s.

4.2. Isotherm: Initial Estimates of b, s. In experimental practices it is common to draw and analyze the isotherm, i.e., the curve of the steady state penetrating flux [bar.J] dependence on the inlet pressure [bar.p], while vacuum pumping is performed at the outlet side. If one targets the Sieverts' law and the (quasi)equilibrium concentration [bar.c] = [GAMMA] [square root of [bar.p]] at the inlet side (hence [bar.J] = D[bar.c]/l), then it is natural to plot the dependence [bar.J] = [bar.J]([square root of [bar.p]]).

Let us analyze the steady state flux balance equation: a [equivalent to] D[(2lb).sup.-1],

[mathematical expression not reproducible]. (25)

Asymptotic analysis shows that the dependence [bar.J] ([square root of [bar.p]]) has a parabolic shape ([bar.J] [varies] [bar.p]) at low inlet pressures [bar.p]:

[mathematical expression not reproducible], (26)

and a straight line form at relatively high pressures [bar.p]:

[a.sup.-2][[GAMMA].sup.2-]p >> 1 [??] [bar.J] = -[D.sup.2][(2[l.sup.2]b).sup.-1] + D[GAMMA][l.sup.-1][square root of [bar.p]]). (27)

Using the straight-line segment of the isotherm we find D[GAMMA]/l (the slope of the straight line), and knowing the estimate of D we determine the initial approximation of the solubility coefficient r. From the intersection of the straight line with the ordinate axis we find b. Knowing the values of [GAMMA] = [square root of 2s[mu]/b] and b we compute s and [PHI] = D[GAMMA].

A graphic illustration (using a minimal required set of the calculated model [[bar.J].sub.1,2,3] values) is presented in Figure 3.

Note that the obtained initial approximations are in good agreement with the original "forgotten" parameters as shown in Table 1.

4.3. The Final Stage of Isothermic Identification. Table 1 reflects only the computational errors arising when solving the direct and inverse problems. The real experimental data are noisy. The identification algorithm uses only integral operators thus ensuring the noise resistance of experimental data treatment. The penetration method is characterized by a significant measurement error, and data on the penetrating flux are required (and this, in turn, requires a more accurate determination of the vacuum system characteristics). The model of the dissolved hydrogen concentration jump at the inlet side is not very precise either. We are brought to a conclusion that the "communicating vessels" stage, where hydrogen pressures are measured over a long time, offers a much higher accuracy of measurements.

Thus, the first stage of the aggregated experiment is perceived as a preliminary estimation of the D, b, s coefficients. It is essential that the solution of the inverse problem of parametric identification is unique, since the results obtained for thin laboratory membranes are extrapolated (recalculated) to the dimensions of real-life structures. The results are "fine-tuned" by means of local variation of the preliminary values of D, b, s in the fast hydrogen permeability model (the ODE system).

5. Thermal Desorption Spectroscopy (TDS)

Pressure values inside the chambers tend to equalize with time. The uniform equilibrium saturation c = [bar.c] = [GAMMA][bar.p], [GAMMA] = [square root of 2[mu]s/b], T = [bar.T] = const corresponds to this pressure p = [bar.p]. We turn off the heating. Next, we create the vacuum at both sides of the membrane and begin slowly reheating the sample from room temperature (the degassing is practically completed) to activate the desorption process.

The thermal desorption flux of hydrogen from the sample is measured by mass spectrometer. What additional information can be obtained from the TDS experiment? The dependence D(T) is assumed to be known (for example, as a result of a series of experiments in the DLR mode at various temperatures). We obtain no new information about s(T) under vacuum conditions since modern, quite powerful vacuum systems allow neglecting the resorption. It remains to more precisely define the two-parameter Arrhenius-like dependence b(T). This is the bulk desorption coefficient (the coefficient of effective recombination of atoms into [H.sub.2] molecules): [J.sub.0,l] = b[c.sup.2.sub.0,l], b = [b.sub.vol].

The heating T(t) is usually linear T(t) = [T.sub.0] + vt. The heating rate v is not too high (<K/s). When the temperature maximum is reached (if degassing is not yet completed), then heating is stopped: T(t) = [T.sub.max]. We refine the boundary conditions taking into account that hydrogen atoms can accumulate at the surface during slow monotone heating and relatively low temperatures.

5.1. Dynamical Boundary Conditions. The surface is a significant potential energy barrier (see [2, p. 177-206]). Hence, the boundary conditions are modeled as follows: c0(t) = c(t, 0), [c.sub.l](t) [equivalent to] c(t,l),

c(0,x) = [bar.c], x [member of] [0,l], t [member of] [0,[t.sub.*]], (28)

[c.sub.0](t) = g (T) [q.sub.0] (t), [c.sub.l](t) = g (T) [q.sub.l] (t), (29)

[[??].sub.0](t) = 2[mu](T)s(T)[p.sub.0](t) - b(T)[q.sup.2.sub.0] (t) + D[[partial derivative].sub.x]c(t,0), (30)

[[??].sub.l](t) = 2[mu](T)s(T)[p.sub.l](t) - b(T)[q.sup.2.sub.l] (t) + D[[partial derivative].sub.x]c(t,l). (31)

Here, [q.sub.0,l](t) are the surface concentrations (x = 0,l, q [equivalent to] dq/dt); g(T) is the parameter of local equilibrium between the surface and the subsurface bulk (coefficient of quick dissolution); b(T) = [b.sub.0] exp{-[E.sub.b]/[RT]} is the surface desorption coefficient:

[J.sub.0,l](t) = b(T(t))[q.sup.2.sub.0,l](t) = b(T(t))[g.sup.-2](T(t))[c.sup.2.sub.0,l](t) [??] b = [b.sub.surface] = [g.sup.2][b.sub.volume]. (32)

The well-posedness of these dynamical boundary conditions was proved in [36].

5.2. "Surface-Bulk" Fluxes Balance. Model (29) of quick dissolution (local equilibrium) on the surface is derived from the more general ratios:

[mathematical expression not reproducible]. (33)

Coefficients [k.sup.-], [k.sup.+] are the descriptors of the rate of dissolution in the bulk and transfer to the surface. When concentrations are nowhere near maximum and D[c.sub.x] [approximately equal to] 0 on the relative scale, we obtain (29), where g = [k.sup.-]/[k.sup.+]. If the surface is isotropic (in terms of [mathematical expression not reproducible]), then the parameter g(T) is a little dependent on temperature. The density of the H atoms adsorption flux can be modeled by the term 2[mu]sp[[1 - [theta]].sup.2] for balance equations (30), (31). For the ranges of weak concentrations and sufficiently high working temperatures the degree of surface coverage satisfies [theta](t) = q(t)/[q.sub.max] << 1. These simplifications are in agreement with the limited information capacities of the TDS experiments. Experimental data are more easily approximated for a large number of parameters. But the uniqueness of the estimations is then failed, and thus essential errors may occur at extrapolation of the results "from thin plate to wall".

5.3. Symmetry Conditions. We shall hereinafter use a contracted notation for simplicity

b(t) [equivalent to] b(T(t)),

D(t) [equivalent to] D(T(t)),

s(t) [equivalent to] s(T(t)). (34)

The generalized quick dissolution coefficient g of local "surface-bulk" equilibrium is assumed to be constant [mathematical expression not reproducible] in the given heating range. The TDS method fulfills the symmetry conditions:

[mathematical expression not reproducible]. (35)

The information capacities of the initial and final stages of the TDS experiment are low. So it is sufficient to limit ourselves to t [less than or equal to] [t.sub.*], when the flux from the sample has decreased by an order of magnitude compared to the maximum. The experimental data are the desorption density curves J(t) or TDS spectra J(T) (T(t) [left and right arrow] t) for different saturation conditions and heating rates. The conversion p(*) [??] J(*) is made more specific by taking into account the parameters of the experimental unit. Modern equipment provides a means of creating deep vacuum ([10.sup.-9]-[10.sup.-7] Torr). For this reason, the component P [equivalent to] 2[mu]sp (P >> 1) is the key control for the saturation stage, but resorption is neglected for the degassing stage (P << 1).

5.4. Diffusion Model with Reversible Capture. We can take into account different channels of diffusion, but the information content of the TDS experiment is limited. Therefore, the coefficient D is assumed to be an integral effective index.

Seeking a write-up of the TDS peaks set, it is handier to use the following model:

[mathematical expression not reproducible], (36)

where [z.sub.v](t, x) is the concentration of hydrogen atoms captured by defects of different types; [a.sup.[+ or -].sub.v] are the coefficients of H capture and release by traps; [Z.sub.v] [equivalent to] [z.sub.v](t,x)/[z.sup.v.sub.max] is the defects saturation degree ([z.sup.v.sub.max] = max [z.sub.v]). Capture is taken into account at its integral level for practical purposes. A more precise definition of the defects' geometry and distribution would add complexity to the model. If, for instance, the defect is not a microcavity but hydride phase inclusions, then at the degassing stage the corresponding coefficient [a.sup.-.sub.j](T) is identically zero and [a.sup.+.sub.j](T) value is positive only if the critical temperature is reached: T(t) [greater than or equal to] [T.sub.crit]. It is easy to simulate the required number of TDS-peaks using different binding energies (coefficients [E.sub.a]). Numerical algorithms based on difference schemes and modeling results were presented in [24, 26]. In this paper we use only (1) (D = [D.sub.eff]). For thin membranes of quickly permeable material used in experiments this approximation is usually sufficient.

5.5. Equilibrium Analysis. Let us take a brief look at the equilibrium ratios at the accepted detail level of modeling. We assume that pressure and temperature are constant. Formally, equilibrium is characterized by all derivatives being equal to zero. Keeping in mind the extensively used Sieverts' law, we shall observe proportionality to the [square root of p] value.

(1) On the surface the following ratio is applied:

2[mu]sp [[1 - [[theta].sub.1]].sup.2] - b[q.sup.2] = 0, [[theta].sub.1] [equivalent to] q/[q.sub.max], [v.sup.2] [equivalent to] 2[mu]s/b[q.sup.2.sub.max]. (37)

Here, [[theta].sub.1] is the degree of surface coverage (next [[theta].sub.i] have a similar meaning). Let us formulate [[theta].sub.1]:

[mathematical expression not reproducible]. (38)

(2) In the "surface-bulk" equilibrium we have ([[theta].sub.2] [equivalent to] c/[c.sub.max], [k.sub.1] [equivalent to] [k.sup.-] [q.sub.max], [k.sub.2] [equivalent to] [k.sup.+] [c.sub.max], [gamma] [equivalent to] [k.sub.1]/[k.sub.2]), then from Formula (33) we obtain

[k.sub.1][[theta].sub.1] [1 - [[theta].sub.2]] = [k.sub.2][[theta].sub.2][1 - [[theta].sub.1] [??] [[theta].sub.2] = [gamma][[theta].sub.1]/1 + [[gamma] - 1] [[theta].sub.1]. (39)

(3) For definiteness, we take into account only one type of traps (m = 1). We use (36). For symmetry, we add the saturation factor [1 - [[theta].sub.2]] for z = [z.sub.v]. Here we have [[theta].sub.3] [equivalent to] z/[z.sub.max], a [equivalent to] [a.sup.+] [z.sub.max]/[[a.sup.-][c.sub.max]]; then we obtain

[[theta].sub.2][1 - [[theta].sub.3] = [[theta].sub.3][1 - [[theta].sub.2] a [??] [[theta].sub.3] = [[theta].sub.2]/a + [1 - a][[theta].sub.2]. (40)

The "saturation-degassing" experiment provides information about a general average concentration [??] in the bulk V = Sl (end surfaces are neglected):

V[??] = 2[[theta].sub.1][q.sub.max]S + [[theta].sub.2][c.sub.max] V + [[theta].sub.3][z.sub.max] V. (41)

Normalize [??] by [c.sub.max] and consider the dependence

[THETA](p) [equivalent to] [??]/[c.sub.max] = 2[q.sub.max]/l[c.sub.max] [[theta].sub.1] + [[theta].sub.2] + [[theta].sub.3] [z.sub.max]/[c.sub.max]. (42)

Let us focus on the curves in the axes ([square root of p], [THETA]). Numerical results are presented in Figure 5 ([bar.a] = [a.sup.+]/[a.sup.-], [c.sub.max] = [10.sup.18], l = 0.1 cm). Quite many parameters have an effect, but one can usually estimate the orders of magnitude to find "the distance to linearity" (to the Sieverts' law adequacy range).

The curves [THETA] = f([square root of p]) have the "growing wave" form. It is noticeable in a wide pressure range only. For the given parameters, the inflection point is the most vivid on the curves marked with pentagrams. The analysis of each additive component in (42) for the total concentration shows that the position of the inflection point is determined by the moment when the first and third addends have turned to the saturation mode with the following prevalent rise of [[theta].sub.2]. Note that the curves for a narrower pressure range are practically linear, in line with ranges of the Sieverts' law. In a wide range of pressures, the wave-like character of the graphs is observed experimentally [37].

5.6. Richardson Approximation. When experimentally estimating the hydrogen permeability of structural materials, Richardson approximation is often used for the penetration flux density:

J(t) = [PHI][[square root of [p.sub.in](t)] - [square root of [p.sub.ou](t)]][l.sup.-1], [p.sub.ou] < [p.sub.in]. (43)

Membranes are usually very thin and their permeability is relatively high. We can accept for the concentration gradient [c.sub.x] = [[c.sub.out] - [c.sub.in]]/l. In accordance with the Fick's law, the formula J = -D[[c.sub.out] - [c.sub.in]]/l is "exact". Here, the function J(t) is the monoatomic hydrogen permeation flux density. The actual boundary volume concentrations cannot however be measured, and thus a quasi-equilibrium approximation is used, where the concentration c is substituted with the equilibrium concentration [bar.c] = [GAMMA][square root of p] in accordance with Sieverts' law (r is the solubility coefficient, or solubility for short). On account of the inequalities [[bar.c].sub.in] > [c.sub.in], [[bar.c].sub.out] < [c.sub.out], this substitution overrates the second factor in the formula J = D[[c.sub.in] - [c.sub.out]]/l : [[bar.c].sub.in] - [[bar.c].sub.out] > [c.sub.in] - [c.sub.out]. Hence the value of D needs to be formally marked down to retain the equality. So, if the value of [PHI] (permeability) is found from experimental data by fitting as indicated in formula (43), then we obtain [PHI] < D[GAMMA]. The [PHI] and [GAMMA] values from the (typically used) formula [PHI] = D[GAMMA] determine the lower bound of the diffusion coefficient D. Furthermore, dissolved diffusing hydrogen is mainly involved in the steady state (quasi-steady state) permeability mode. In the context of the model (29)-(31), we obtain [bar.c] = [GAMMA] [square root of p], [GAMMA] = g [square root of 2[mu]s/b]. The "saturation-degassing" experiment yields the total value [??] > [bar.c] and an overstated (for the permeability problem) value [[GAMMA].sub.max] : [??] = [[GAMMA].sub.max] [square root of p].

Hence, we can estimate the [PHI] value by fitting using formula (43). This information has practical value as a convenient coefficient for translation from pressures to flux. If, however, the D value is taken from one experiment (or paper) and the [GAMMA] value is taken from another source, then, strictly speaking, we get a ranking of three different numbers [PHI] < D[GAMMA] < D[[GAMMA].sub.max]. If the material has a high level of trapping by defects, then the calculated permeability can be an order of magnitude higher than the real permeability. The permeability coefficient [PHI] (as a parameter of (43)) has an S-form (Arrhenius-like) of the saturation curve based on the order of pressures. It is only for relatively high pressures (where [c.sub.0,l] are near Sieverts' concentrations) that we get [PHI] [approximately equal to] D[GAMMA].

6. Functional Differential TDS Equation

Identification of TDS spectra is required not only to reveal the causes of different thermal desorption peaks, but also to enable numerical extrapolation and generalization of the results received for laboratory samples (l usually is fractions of mm). Model (36) gives the possibility to get any number of peaks using traps with different parameters [a.sup.[+ or -].sub.v]. The question, however, is whether different peaks can occur when degassing an almost homogeneous material. To answer this question let us restrict ourselves to the basic diffusion equation (1), but retaining symmetric dynamical boundary conditions (29)-(31) and (35). The surface is considered isotropic in terms of g = const over the heating range. The resorption during vacuum building is neglected. Thus, we are limited to a minimal number of parameters for the model which takes into account the dynamical interplay of surface processes and diffusion. In the following, let us operate at this approximation.

The comparison of simulated and experimental TDS spectra with a focus on parametric identification requires only the surface concentration (J = b[q.sup.2]). It is reasonable to try to avoid iterative solution of the boundary value problem for interim approximations of the model parameters [D.sub.0], [E.sub.D], [b.sub.0], [E.sub.b], [s.sub.0], [E.sub.s], g. To this end, we will run the transformations to reduce the problem to the integration of a low order ODE system.

6.1. Derivation of Riccati-Type Equation. The accepted TDS degassing model is [[partial derivative].sub.t]c = D(t)[[[partial derivative].sup.2.sub.x]c, c(0, x) = [bar.c],

[mathematical expression not reproducible]. (44)

Let us replace the time t' = [[integral].sup.t.sub.0] D d[tau] (keeping the former notation t):

[[partial derivative].sub.t]c (t,x) = [[partial derivative].sup.2.sub.x]c (t,x), c(0,x) = [bar.c], (45)

[c.sub.0,l] = gq(t), [[partial derivative].sub.x][c|.sub.0] = - [[partial derivative].sub.x][c}.sub.l] = [??](t) + J(t)[D.sup.-1](t). (46)

Here q(t) is considered as the parameter and (46) is an additional relation for the linear problem (45). We perform a replacement to get homogenous boundary conditions:

[mathematical expression not reproducible]. (47)

Let us write the solution of the linear boundaryvalue problem using the source function (Green's function):

[mathematical expression not reproducible]. (48)

Boundary conditions contain the derivative [[??].sub.x](t, 0):

[mathematical expression not reproducible]. (49)

For [tau] = t we have divergent series. So, term-by-term integration is implied. For the original time t we get [[partial derivative].sub.x]c(t, 0) = [[partial derivative].sub.x][??](t,0) = [[partial derivative].sub.x]c(t,l) = [[partial derivative].sub.x][??](t,l),

[[partial derivative].sub.x]c(t,0)

= - 4g/l [summation]' [[integral].sup.t.sub.0] [??]([tau]) exp {[n.sup.2][[pi].sup.2]/[l.sup.2] [[integral].sup.t.sub.[tau]] D(s) ds} d[tau]. (50)

Finally, the dynamic boundary condition is written in the integrodifferential form:

[mathematical expression not reproducible]. (51)

The resultant equation (51) with quadratic nonlinearity will be classified as a functional differential Riccati equation of the neutral type. The equation is equivalent to the original boundary value problem in that the solution q(t) uniquely determines the solution c(t, x). The analogy with the functional differential equation of the form x(t) = F[t, x(t), [??](t), [x.sub.[0,t]],[[??].sub.[0,t]] of the neutral type [38] is that it is impossible to eliminate the derivative q on the right side of the equation through integration by parts lest a divergent series arises. We are concerned with the time interval [[t.sub.1],[t.sub.2]] [subset] (0,[t.sub.*]), which corresponds to the TDS peak. Measurement for t [approximately equal to] 0,[t.sub.*], yields little information. There is a voluminous body of literature on Riccati-type equations (including matrix equations for the optimal control theory).

6.2. Dimensionless Form of the Problem. For more comfortable modeling we turn to dimensionless variables using substitution rules: t' = [[integral].sup.t.sub.0] D([tau])d[tau]/[l.sup.2], x' = x/l, v = q/[bar.q] ([bar.c] = g[bar.q]). Retaining the notation t, we obtain

[mathematical expression not reproducible], (52)

where [??](t) [equivalent to] [bar.q]b(t)[l.sup.2]/D(t) is a dimensionless parameter of quadratic "desorption".

6.3. Initial Data. We specify the factor [??](t) for quadratic nonlinearity. Initial saturation is conducted under relatively high temperature T = [bar.T] = const and pressure [bar.p] = const. After the steady state of saturation is attained we get

[mathematical expression not reproducible]. (53)

Hereafter, let us be guided by a maximum limit of surface concentration at around [q.sub.max] ~ [10.sup.15]-[10.sup.16] (monolayer in the context of geometric statics). If during initial saturation the surface concentration for the given model parameters is [bar.q] > [10.sup.14], the degree of surface coverage must be taken into account. Then, to calculate the initial value of [bar.q] we take the ratio 2[mu]s[bar.p][[1 - [bar.q][q.sup.-1.sub.max]].sup.2] = b[[bar.q].sup.2] instead of 2[mu]s[bar.p] = b[[bar.q].sup.2].

In the meantime, this a priori restriction (arising from the assumption that stationary "balls" are ordered geometrically on a plane) is highly questionable. For a dynamics model, it appears that the concentration threshold [q.sub.max] may be higher, if it is specified what meaning is attached to the term "bubbling" surface layer (at a quite high temperature). The real [bar.q] value is somewhat conventional, since it is strongly influenced by the experiment pretreatment (building of vacuum before the start of heating). However, most of the hydrogen is in the bulk, the diffusion equation (and the diffusion process itself) has a smoothing effect, and we are interested in evident thermal desorption peaks, since the initial and final time intervals of the experiment offer little information. For initial calculations it is therefore sufficient to correctly estimate the magnitude of the "effective" [bar.q] concentration. High precision requirements are noncritical here.

7. Extracting Integrable Singularity

The functional differential equation of thermal desorption (52) has a singularity hindering its numerical solution. The function,

[THETA](s) = 4 [summation]' exp{-[n.sup.2][[pi].sup.2]s}, [summation]' [equivalent to] [summation over (n=1,3,5,...)], (54)

takes finite values for s >0. The series is rapidly converging for large s. Formally, substituting s = 0 (the integration variable t reaches the upper limit t) we obtain a divergent series. This can be "fixed" by term-by-term integration. The order of terms in the series is O([n.sup.-2]) (n = 2i - 1, i [greater than or equal to] 1), wherefore convergence is slow. Let us formulate the problem. The equation has to be approximated by a low order ODE system to enable application of standard software packages (for example, Scilab).

Let us run some transformations using the theory of Jacobi theta functions to explicitly extract the integrable singularity. We consider the presentation (14):

[[theta].sub.3](t,0) = 1 + 2 [[infinity].summation over (n=1)] exp{-[n.sup.2][[pi].sup.2]t} =1/[square root of [pi]t] [[infinity].summation over (-[infinity])] exp {-[n.sup.2]/t}. (55)

The series on the left is rapidly converging for large t. But the series on the right is rapidly converging for small t. Let us run auxiliary transformations:

[mathematical expression not reproducible]. (56)

Proceeding from here, the last series is subtracted from the first series and the result is doubled. The following expression is obtained for s > 0:

[mathematical expression not reproducible], (57)

q [equivalent to] - exp{-1/(4s)}. The series Q is very rapidly converging for small s. As s [right arrow] +0 we have Q [right arrow] 0 and the integrable singularity [THETA] [approximately equal to] 1/[square root of [pi]s]. The graph for Q(s) has an S-shaped ("Arrhenius-like") saturation curve form and S(1) [approximately equal to] 0.9996. The function Q/[square root of [pi]s] increases monotonically to a maximum ([approximately equal to] 0.828 for s [approximately equal to] 0.334) and then decreases monotonically. To represent the series Q(s) it is reasonable to use only a small number (5-8) of series terms. As an alternative to (52), using the above-described presentation of [THETA](s) we obtain the equation

[mathematical expression not reproducible], (58)

where the fraction (weak singularity under the integral) rapidly decreases from infinity ([tau] = t) to near zero ([tau] = t-1). For t > 1, the lower limit of the integral can be replaced with t-1. Thus, the neglectable background can be easily identified in the original physical time using the relation [mathematical expression not reproducible]. The compact functional differential ([??] [equivalent to] dv/dt) TDS equation (58) with initial data v(0) = 1 replaces the original nonlinear boundary value problem (from Section 6.1) with dynamic boundary conditions in the sense that formally only the dynamics of the surface concentration (desorption) is required for TDS spectrum construction. Note that equation (58) contains a Caputo fractional derivative [d.sup.1/2] v(t)/d[t.sup.1/2] = [[integral].sup.t.sub.0][??]([tau])[[pi](t - [tau])].sup.-1/2] d[tau]. The theory of integrodifferential equations (the "residual" integral containing the function Q can be transformed by parts) is a surging field of modern mathematics and its applications.

8. Numerical Method and Computer Simulation

To be specific in the paper, we use data for nickel and steel (12Cr18Ni10Ti) [6], tungsten [5], and beryllium [39,40]. Estimates depend substantially on the experimental conditions and sample pretreatment. So the values are perceived as a model for numerical illustrations. The parameters common for all the materials are l = 0.1 cm, T0 = 300 K, [E] = J/mol. The assumed values of model parameters are the following: (steel) [b.sub.0] = 1.28 x [10.sup.-9] [cm.sup.2]/s, [E.sub.b] = 80 x [10.sup.3], [D.sub.0] = 3.09 x [10.sup.-4] [cm.sup.2]/s, [E.sub.D] = 31 x [10.sup.3], g = 50 [cm.sup.-1], [E.sub.g] = -5 x [10.sup.3], [s.sub.0] = 0.6, [E.sub.s] = 110 x [10.sup.3], [bar.c] = 6 x [10.sup.17] at.H/[cm.sup.3], [bar.T] = 1170 K, [??] = 1 K/s; (Ni) [b.sub.0] = 1.53 x [10.sup.-14], [E.sub.b] = 43.2 x [10.sup.3], [D.sub.0] = 7.5 x [10.sup.-3], [E.sub.D] = 40 x [10.sup.3], g = 100, [E.sub.g] = 0, [s.sub.0] = 1.8 x [10.sup.-2], [E.sub.s] = 61.4 x [10.sup.3], [bar.p] = 37.4 Torr, [bar.T] = 770, [??] = 0.5; (W) [b.sub.0] = 6 x [10.sup.-4], [E.sub.b] = 69.559 x [10.sup.3], [D.sub.0] = 4.1 x [10.sup.-3], [E.sub.D] = 37.629 x [10.sup.3], [g.sub.0] = [10.sup.4], [E.sub.g] = 15 x [10.sup.3], [s.sub.0] = 9 x [10.sup.-3], [E.sub.s] = 15, [bar.p] = 670, [bar.T] = 1300, [??] =5; (Be) [b.sub.0] = 3.08 x [10.sup.-9], [E.sub.b] = 57.43 x [10.sup.3], [D.sub.0] = 3 x [10.sup.-3], [E.sub.D] = 28 x [10.sup.3], [g.sub.0] = 200, [E.sub.g] = 1.5 x [10.sup.3], [s.sub.0] = 1.44 x [10.sup.-4], [E.sub.s] = 1.82 x [10.sup.3], [bar.p] = 760, [bar.T] = 1150, [??] = 5.

The main role in the degassing dynamics belongs to quadratic desorption. We therefore approximate the integral term in (58). The after-effect horizon here is h < 1 ([tau] [member of] [t, t-h] for dimensionless time, Q(1) [approximately equal to] 0.999). Let us fix h~ 0.3-0.4 due to the smooth function Q(s)/[square root of s] graph (Figure 4) and the trapezoid rule for numerical integration. The replacement of t' = [[integral].sup.t.sub.0] Dd[tau]/[l.sup.2] is naturally targeted at the diffusion time scale [l.sup.2]/D so that the step h corresponds to a significant segment of the experiment physical time. Then, TDS equation (58) can be approximated step-by-step by a low order ODE system on segments of dimensionless time of length h.

The greatest contribution to the integral is made by the value [??]([tau]), [tau] [approximately equal to] t due to nonlimited (but integrable) singularity. Thus, the quadratic approximation

[??](t) [approximately equal to] [??](t) + A[t - [tau]] + B[[t - [tau]].sup.2] (59)

is used due to function V concavity. Let us consider the current segment of time t [member of] [kh, (k + 1)h], k [greater than or equal to] 0. The conditions

[mathematical expression not reproducible] (60)

determine the values of A(t), B(t) (constants on [tau]):

[mathematical expression not reproducible], (61)

[tau] [member of] [kh,t] (t [member of] [kh,(k + 1)h]). Represent the integral from Formula (58) as a sum [tau] [member of] [0, kh], [tau] [member of] [kh, t]. For the second additive, the singularity [??](t)/[square root of t - [tau]] is explicitly integrated by substituting (9). The trapezoid rule and mean value theorem are used for the integral without singularity (Q(+0)/[square root of +0] = 0)

[mathematical expression not reproducible]. (62)

To approximate the first integral

[mathematical expression not reproducible] (63)

we use only several terms (for definiteness n = 1,3) on the right-hand side taking into account the factor [n.sup.2]. The negative additives ([??] < 0) thereby drop out. Let us compensate for that by replacing kh with t and introducing the additional variables [w.sub.1,2](t):

[mathematical expression not reproducible]. (64)

As a result we obtain an ODE system instead of (58):

[mathematical expression not reproducible]. (65)

Ellipsis stands for the right-hand side of the first line of the first equation (-[??] ... [??](kh)). The series for Q(s) converges very rapidly. It is sufficient to use 5-8 terms for software.

The sought function v(t) = q(t)/[bar.q] on the current segment of dimensionless time (in the system (65) t = t', t' [member of] [kh, (k + 1)h], k [greater than or equal to] 0) is computed by the Runge-Kutta 4th-order method for integration of ODE systems (the authors used Scilab software). If the TDS spectrum has two peaks, then [w.sub.3,4] should be additionally used to improve the modeling accuracy. Returning to the physical time, we get the model spectrum J(T) = b(T)[q.sup.2](t) (t' [left and right arrow] t [left and right arrow] T(t)).

Qualitatively, thermal desorption spectra of metal structural materials have a typical form. Figure 6 illustrates numerically simulated TDS spectra for the above-listed hydrogen permeability parameter values of the structural materials at different heating rates. Even when defects (traps) are not taken into account, the suggested model can yield different types of two-peak graphs (where the low temperature peak is more pronounced or where the peaks are comparable). For reference, the change of the so-called transport parameter W = l[b.sub.vol][bar.c]/D = lb[bar.c](D[g.sup.2]) is shown in the top part of the graphs. This parameter is crucial for the study of membrane hydrogen permeability [5]. We supposed that SLR corresponds to W < 10-2 and DLR corresponds to W > [10.sup.4]. For beryllium and steel, the low temperature TDS peak takes place in the range where diffusion and surface processes play commensurable roles. The high temperature TDS peak occurs in DLR. For nickel, surface processes and diffusion have similar effect throughout the experimental temperature range, and only a peak-like step appears in the TDS spectrum at a low temperature.

Let us focus on TDS spectra for tungsten because these spectra are qualitatively different from the previous ones. The aforementioned parameter values "for tungsten" are, of course, formal. Besides microimpurities, the parameters depend on how the sample surface was treated. This is especially true because in practice a plate is thin, the volume is small, and the effect of surface processes is more vivid. This is one of the reasons for such a high variation of the quantitative estimates of hydrogen permeability parameters. Another reason (but not the last) is the following. Different models are used for experimental data treatment (although coefficients formally have the same name). For the model data accepted here, a narrow splash followed by a lengthy movement to the second peak (less visible) was observed. In this model the sample has no high-capacity traps (standard diffusion equation), but more detailed consideration is given to the surface (see (29)-(31), the plate is thin). At first, near-surface hydrogen is actively desorbed. Then, diffusion is slowly activated by heating. The concentration gradient is substantial, and pumping to the surface is growing. These effects are usually categorically explained by the presence of traps, not by the dynamics of "surface-bulk" interplay. An amazingly similar experimental curve (marked with black circles) is found in [41], Figure 1, although this paper discussed deuterium implantation on tungsten. The essence of our model results is that various TDS peaks may have other causes apart from the routinely blamed traps.

Figure 7 separately shows the dynamics of surface processes and diffusion. High-temperature peaks correspond to the significantly enlarged desorption coefficient and activation of the diffusion afflux from the bulk under heating. Peaks at a low temperature take place where the diffusion towards the surface is not high but near-surface concentrations are higher.

Let us briefly present the numerical modeling algorithm.

(i) Set the parameters [s.sub.0], [E.sub.s], [b.sub.0], [E.sub.b], [D.sub.0], [E.sub.D], g. Determine the values of [bar.c], [??](t'), n = gl under the saturation conditions [bar.p], [bar.T] as described in Initial saturation stage. If the equilibrium bulk concentration [bar.c] of dissolved diffusing hydrogen is known, then presetting of s([bar.T]) is not needed. Other scenarios of initial saturation are possible (see discussion in Initial saturation stage). The required adjustment of initial data is determined by the specifics of the actual TDS experiment.

(ii) The low order ODE system of type (65) is numerically integrated step by step in dimensionless time (on an h length interval). Other ODE approximations were presented in [27]. This procedure is standard for every modern software.

(iii) The graph J(t) = b(T(t))[q.sup.2](t) or spectrum J(T) is plotted in physical time using the function v(t') = q(t')/[bar.q].

9. Inverse Problem of Parametric Identification

The presented algorithm of numerical modeling allows quick scanning of different scenarios and operating conditions of a material (including the heating law and extrapolation of the results with l increase). This statistical information is useful when designing the strategy of experimental research. New materials (various alloys) first have to be analyzed for their hydrogen permeability. In dealing with this task we encounter inverse mathematical physics problems, which are well-known for their difficulty. Let us suppose that the diffusion coefficient D(T) is known (Daines-Berrer method is usually used for the DLR of permeability). Let us formulate the problem: to estimate the surface parameters [b.sub.0], [E.sub.b], g using the J(t) (desorption flux) information. These conditions correspond to the real conditions of an experiment where desorption dynamics closely correlates with diffusion in the bulk. The problem of estimating the desorption coefficient ([b.sub.vol] = b/[g.sup.2]) for the situation where accumulation on the surface can be neglected (q [approximately equal to] 0) was presented in [25].

Since the function D(t) [equivalent to] D(T(t)) is known, it is reasonable to move directly to dimensionless time t', which is oriented at the characteristic time of the diffusion transfer [l.sup.2]/D. The old notation t is retained not to complicate the formulas. Let us present the nonlinear additive term in the system (65) in the following form:

[mathematical expression not reproducible]. (66)

Here, [alpha] is a function of the parameter k. It is obtained through transformations using notations from (52) and (65). Elementary but somewhat lengthy formulas are omitted. The function J(t) (t = t') is known from measurements. After substitution into (65) we obtain a system of three linear ODE. It is reasonable to add the variables [w.sub.3,4] to improve the calculation accuracy. The measurements are usually noisy, but the function J(t) comprised of the right-hand side of the ODE and is smoothed by integration.

Note that the right-hand sides of the equations now depend only on one estimated dissolution parameter k = gl. The identification algorithm can thus be decomposed. By setting the value of g we can eventually calculate the dependence q(t) = q(t; g) in the original physical time t. Be reminded now that

[mathematical expression not reproducible], (67)

where [psi](t) [equivalent to] J/[q.sup.2], T = T(t). Since heating is monotonic t [left and right arrow] T(t), we can introduce the coordinates X [equivalent to] [10.sup.3]/T, Y [equivalent to] ln [psi]. The parametric curve X(t), Y(t) is obtained on the plane (X, Y). Judging by the ratio (67), this curve has to be a straight line with a negative slope. Hence, we have the criterion for the choice of the "correct" value for g: the curve X(t), Y(t) = ln{J(t)/[q.sup.2](t; g)} has to be a straight line segment on the plane (X, Y). This subproblem is scalar, where only g varies.

Formally, we prolong the straight line segment until it crosses the coordinate axes. The crossing with the y-axis (X = 0) yields the ln [b.sub.0] value. The crossing with the abscissa [bar.x] = [10.sup.3] x R ln [b.sub.0]/[E.sub.b] determines the [E.sub.b] value.

The algorithm is laborious due to iterative use of the procedure of numerical solving of the initial problem for the ODE system. It requires some effort and familiarity with mathematical packages, but it is much easier and quicker to use the standard built-in operation than to perform iterative solution of the original nonlinear boundary value problem with three varying parameters.

The presented linearization method is illustrated in Figure 8. Model spectra were plotted to test the algorithm. A system of five equations (for the variables v(t), [w.sub.1-4](t)) was numerically solved to obtain more accurate estimates for two-peak spectra. In addition, a series of experiments was performed for spectra with a random error not higher than 20% (we used standard Scilab uniform random number generator). The identification algorithm based on ODE system integration demonstrated the noise resistance of experimental data treatment.

The initial value of g yields the mass balance

[mathematical expression not reproducible]. (68)

The experiments have demonstrated that the accuracy of such estimation decreases for g [approximately equal to] 1. Curves based on the "correct" numerical k value (the resultant curve is close to the straight line segment) and curves for 20% deviation from k are shown. The error of the identification algorithm testing is less than several percentages. This error appears due to the numerical error of direct and inverse problem solving. There is high sensitivity (appearing as vertical "beak" singularity) to situations where the "true" k value is exceeded. Mathematical singularity appears because, formally, the function v(t) changes sign, taking also negative values. We did not take efforts to avoid this "nonphysical" transition because it is a vivid sign of a wrong k value. The same numerical experiments were done for noisy spectra. Activation energies are determined more accurately because energy parameters more actively influence desorption during heating. From the mathematical point of view, there is no need to strive for visual agreement of simulated and experimental curves (because the experimental error is tens of percentages). ODE approximation is quite adequate.

10. Conclusions

It is advisable to aggregate the thermal desorption (TDS) and penetration experiments (with and without vacuum pumping) to make the measurements substantially more informative for further estimation of the hydrogen permeability parameters and to improve the accuracy of parameter identification. This paper suggests a cascade experiment technique and the corresponding mathematical software.

The identification algorithm uses only integral operators thus ensuring the noise resistance of experimental data treatment. The penetration method with vacuum pumping is characterized by a significant measurement error, and data on the penetrating flux are required (and this, in turn, requires a more accurate determination of the vacuum system characteristics). The model of the dissolved hydrogen concentration jump at the inlet side is not very precise either. We are brought to a conclusion that the "communicating vessels" stage, where molecular hydrogen pressures are measured over a long time, is characterized by a much higher accuracy of measurements.

The first stage of the aggregated experiment is perceived as preliminary estimation of the diffusion D, desorption b = [b.sub.vol], and absorption s coefficients. It is essential that the solution of the inverse problem of parametric identification is unique, since the results obtained for thin laboratory membranes are extrapolated (recalculated) to the dimensions of real-life structures. The results are "fine-tuned" by means of local variation of the preliminary values of D, b, s in the ODE-model of fast hydrogen permeability.

The final stage deals with the problem of identification of the spectra of hydrogen thermal desorption. This problem is of high relevance for the nuclear power industry. Qualitatively, the identification consists of revealing the causes of desorption peaks. It is usually assumed that TDS peaks appear due to hydrogen release from traps with different binding energies. Here, it was demonstrated using a diffusion model for a homogeneous material that if surface processes are taken into account, two-peak spectra can be obtained even for very thin experimental samples. The tendency to resort to the "theory of different traps" as the only explanation is understandable, but the volume of our samples was near zero for the trapping capacity to manifest itself.

The nonlinear boundary value problem (standard diffusion equation with dynamical boundary conditions) is reduced to the functional differential equation for the surface concentration, because nothing but desorption dynamics is required to plot a TDS spectrum. An effective numerical algorithm oriented to the use of mathematical packages (including freeware) is proposed. The main final output of this part of the paper is a geometrically transparent method for solving the inverse problem of surface parameter identification where desorption and diffusion in the bulk are dynamically interrelated.

https://doi.org/10.1155/2018/4628346

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this article.

Authors' Contributions

All authors contributed equally to the writing of this article. All authors read and approved the final manuscript.

Acknowledgments

The study was carried out under state order to the Karelian Research Centre of the Russian Academy of Sciences (Institute of Applied Mathematical Research KarRC RAS).

References

[1] G. Alefeld and J. Volk, Eds., Hydrogen in Metals, Springer-Verlag, Berlin, Germany, 1978.

[2] A. P. Zakharov, Ed., Interactions of Hydrogen with Metals, Nauka, Moscow, Russia, 1987.

[3] M. Ball and M. Wietschel, Eds., The Hydrogen Economy: Opportunities and Challenges, Cambridge University Press, New York, NY, USA, 2009.

[4] M. Hirscher, Ed., Handbook of Hydrogen Storage: New Materials for Future Energy Storage, Wiley-VCH, Weinheim, Germany, 2010.

[5] A. A. Pisarev, I. V. Tsvetkov, E. D. Marenkov, and S. S. Yarko, Hydrogen Permeability through Metals, MEPhI, Moscow, Russia, 2008.

[6] A. A. Yukhimchuk, Ed., Interaction of Hydrogen Isotopes with Structural Materials, Russian Federal Nuclear Center, Sarov, Russia, 2008.

[7] R. A. Varin, T. Czujko, and Z. S. Wronski, Nanomaterials for Solid State Hydrogen Storage, Springer, New York, NY, USA, 2009.

[8] M. V. Lototskyy, V. A. Yartys, B. G. Pollet, and R. C. Bowman Jr., "Metal hydride hydrogen compressors: A review," International Journal of Hydrogen Energy, vol. 39, no. 11, pp. 5818-5851, 2014.

[9] Y. P. Cherdantsev, I. P. Chernov, and Y. I. Tyurin, Methods of Studying Metal-Hydrogen Systems, TPU, Tomsk, Russia, 2008.

[10] I. E. Gabis, "The method of concentration pulses for studying hydrogen transport in solids," Technical Physics, vol. 44, no. 1, pp. 90-94, 1999.

[11] M. D. Dolan, "Non-Pd BCC alloy membranes for industrial hydrogen separation," Journal of Membrane Science, vol. 362, no. 1-2, pp. 12-28, 2010.

[12] S. Kozhakhmetov, N. Sidorov, V. Piven, I. Sipatov, I. Gabis, and B. Arinov, "Alloys based on Group 5 metals for hydrogen purification membranes," Journal of Alloys and Compounds, vol. 645, no. 1, Article ID 33306, pp. S36-S40, 2015.

[13] A. Voyt, N. Sidorov, I. Sipatov, M. Dobrotvorskii, V Piven, and I. Gabis, "Hydrogen solubility in V85Ni15 alloy," International Journal of Hydrogen Energy, vol. 42, no. 5, pp. 3058-3063, 2017.

[14] F. J. Castro and G. Meyer, "Thermal desorption spectroscopy (TDS) method for hydrogen desorption characterization (I): theoretical aspects," Journal of Alloys and Compounds, vol. 330-332, pp. 59-63, 2002.

[15] D. A. Indeitsev and B. N. Semenov, "About a model of structural-phase transformations under hydrogen influence," Acta Mechanica, vol. 195, no. 1-4, pp. 295-304, 2008.

[16] E. Evard, I. Gabis, and V. A. Yartys, "Kinetics of hydrogen evolution from MgH2: Experimental studies, mechanism and modelling," International Journal of Hydrogen Energy, vol. 35, no. 17, pp. 9060-9069, 2010.

[17] Y. V. Zaika and N. I. Rodchenkova, "Boundary-value problem with moving bounds and dynamic boundary conditions: diffusion peak of TDS-spectrum of dehydriding," Applied Mathematical Modelling: Simulation and Computation for Engineering and Environmental Systems, vol. 33, no. 10, pp. 3776-3791, 2009.

[18] K. Schmid, V. Rieger, and A. Manhard, "Comparison of hydrogen retention in W and W/Ta alloys," Journal of Nuclear Materials, vol. 426, no. 1-3, pp. 247-253, 2012.

[19] E. Legrand, A. Oudriss, C. Savall, J. Bouhattate, and X. Feaugas, "Towards a better understanding of hydrogen measurements obtained by thermal desorption spectroscopy using FEM modeling," International Journal of Hydrogen Energy, vol. 40, no. 6, pp. 2871-2881, 2015.

[20] A. K. Belyaev, A. M. Polyanskiy, V. A. Polyanskiy, C. Sommitsch, and Y. A. Yakovlev, "Multichannel diffusion vs TDS model on example of energy spectra of bound hydrogen in 34CrNiMo6 steel after a typical heat treatment," International Journal of Hydrogen Energy, vol. 41, no. 20, pp. 8627-8634, 2016.

[21] E. A. Denisov, M. V. Kompaniets, T. N. Kompaniets, and I. S. Bobkova, "Peculiarities of hydrogen permeation through Zr-1%Nb alloy and evaluation of terminal solid solubility," Journal of Nuclear Materials, vol. 472, pp. 13-19, 2016.

[22] Y. V. Zaika and E. P. Bormatova, "Parametric identification of hydrogen permeability model by delay times and conjugate equations," International Journal of Hydrogen Energy, vol. 36, no. 1, pp. 1295-1305, 2011.

[23] Yu. V. Zaika and N. I. Rodchenkova, "Hydrogen-solid boundary-value problems with dynamical conditions on surface," in Mathematical Modelling, C. R. Brennan, Ed., pp. 269-302, Nova Science Publishers, New York, NY, USA, 2013.

[24] N. I. Rodchenkova and Y. V. Zaika, "Numerical modelling of hydrogen desorption from cylindrical surface," International Journal of Hydrogen Energy, vol. 36, no. 1, pp. 1239-1247, 2011.

[25] Y. V. Zaika and E. K. Kostikova, "Determination of effective recombination coefficient by thermodesorption method," International Journal of Hydrogen Energy, vol. 39, no. 28, pp. 15819-15826, 2014.

[26] Y. V. Zaika and E. K. Kostikova, "Computer simulation of hydrogen thermodesorption," Advances in Materials Science and Applications, vol. 3, no. 3, pp. 120-129, 2014.

[27] Y. V. Zaika and E. K. Kostikova, "Computer simulation of hydrogen thermal desorption by ODE-approximation," International Journal of Hydrogen Energy, vol. 42, no. 1, pp. 405-415, 2017.

[28] D. Y. Andronov, D. G. Arseniev, A. M. Polyanskiy, V. A. Polyanskiy, and Y. A. Yakovlev, "Application of multichannel diffusion model to analysis of hydrogen measurements in solid," International Journal of Hydrogen Energy, vol. 42, no. 1, pp. 699-710, 2017.

[29] G. Song, M. D. Dolan, M. E. Kellam, D. Liang, and S. Zambelli, "V-Ni-Ti multi-phase alloy membranes for hydrogen purification," Journal of Alloys and Compounds, vol. 509, no. 38, pp. 9322-9328, 2011.

[30] K. A. Terrani, M. Balooch, D. Wongsawaeng, S. Jaiyen, and D. R. Olander, "The kinetics of hydrogen desorption from and adsorption on zirconium hydride," Journal of Nuclear Materials, vol. 397, no. 1-3, pp. 61-68, 2010.

[31] Y. Zhang, R. Maeda, M. Komaki, and C. Nishimura, "Hydrogen permeation and diffusion of metallic composite membranes," Journal of Membrane Science, vol. 269, no. 1-2, pp. 60-65, 2006.

[32] V. N. Alimov, A. O. Busnyuk, M. E. Notkin, and A. I. Livshits, "Hydrogen transport by group 5 metals: Achieving the maximal flux density through a vanadium membrane," Technical Physics Letters, vol. 40, no. 3, pp. 228-230, 2014.

[33] Yu. V. Zaika, N. I. Rodchenkova, and N. I. Sidorov, "Modeling of H2-permeability of alloys for gas separation membranes," Computer Research and Modeling, vol. 8, no. 1, pp. 121-135, 2016 (Russian).

[34] E. T. Whittaker and G. N. Watson, A Course of Modern Analysis, Cambridge University Press, New York, NY, USA, 1996.

[35] S. Lang, Elliptic Functions, Addison-Wesley Educational Publishers, 1974.

[36] Y. V. Zaika, "Solvability of the equations of a model of gas transport across a membrane with dynamic boundary conditions," Computer Mathematics and Mathematical Physics, vol. 36, no. 12, pp. 1731-1741, 1996.

[37] H. Yukawa, G. X. Zhang, N. Watanabe, M. Morinaga, T. Nambu, and Y. Matsumoto, "Analysis of hydrogen diffusion coefficient during hydrogen permeation through niobium and its alloys," Journal of Alloys and Compounds, vol. 476, no. 1-2, pp. 102-106, 2009.

[38] J. K. Hale, Theory of Functional Differential Equations, Springer, New York, NY, USA, 1977.

[39] A. V. Samsonov, A. Y. Koren'kov, I. E. Gabis, and A. A. Kurdyumov, "Limiting role of desorption in hydrogen transport across a deposited beryllium film," Technical Physics, vol. 43, no. 1, pp. 114-116, 1998.

[40] I. L. Tazhibaeva, V. P. Shestakov, A. K. Klepikovet al., "Modeling of hydrogen release from irradiated beryllium," Bulletin of National Nuclear Center of the Republic of Kazakhstan, vol. 1, pp. 37-41, 2000 (Russian).

[41] Y. Oya, Y. Inagaki, S. Suzuki et al., "Behavior of hydrogen isotope retention in carbon implanted tungsten," Journal of Nuclear Materials, vol. 390-391, no. 1, pp. 622-625, 2009.

Yury V. Zaika (iD), Natalia I. Rodchenkova (iD) and Ekaterina K. Kostikova (iD)

Institute of Applied Mathematical Research of the Karelian Research Centre of the Russian Academy of Sciences, 11 Pushkinskaya St., Petrozavodsk 185910, Russia

Correspondence should be addressed to Yury V. Zaika; zaika@krc.karelia.ru

Received 7 December 2017; Accepted 2 May 2018; Published 4 June 2018

Academic Editor: Alexander Iomin

Caption: Figure 1: The establishment of fluxes J(t).

Caption: Figure 2: Dynamics of pressures [p.sub.0,l](t).

Caption: Figure 3: Extrapolation of the isotherm of [bar.J].

Caption: Figure 4: Functions Q(s) and Q(s)/[square root of [pi]s].

Caption: Figure 5: Saturation degree.

Caption: Figure 6: TDS-spectra.

Caption: Figure 7: Dynamics of surface processes and diffusion.

Caption: Figure 8: Estimation of desorption and dissolution parameters.

Table 1 Original value Approximate value D 2.00 x [10.sup.-5] 2.0012 x [10.sup.-5] b 5.72 x [10.sup.-24] 5.8371 x [10.sup.-24] s 1.20 x [10.sup.-4] 1.2221 x [10.sup.-4] [GAMMA] 2.00 x [10.sup.20] 1.9984 x [10.sup.20] [PHI] 4.00 x [10.sup.15] 3.9993 x [10.sup.15]

Printer friendly Cite/link Email Feedback | |

Title Annotation: | Research Article |
---|---|

Author: | Zaika, Yury V.; Rodchenkova, Natalia I.; Kostikova, Ekaterina K. |

Publication: | Advances in Mathematical Physics |

Article Type: | Report |

Geographic Code: | 1USA |

Date: | Jan 1, 2018 |

Words: | 13893 |

Previous Article: | Heat Transport and Majorana Fermions in a Superconducting Dot-Wire System: An Exact Solution. |

Next Article: | Exact Traveling Wave Solutions of Certain Nonlinear Partial Differential Equations Using the (G'/[G.sup.2])-Expansion Method. |

Topics: |