Abstract

One of the technological challenges for hydrogen materials science is the currently active search for structural materials with important applications (including the ITER project and gas-separation plants). 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). The article presents boundary value problems of hydrogen permeability and thermal desorption with dynamical boundary conditions. A numerical method is developed for TDS spectrum simulation, where only integration of a nonlinear system of low order ordinary differential equations is required. The main final output of the article is a noise-resistant algorithm for solving the inverse problem of parametric identification for the aggregated experiment where desorption and diffusion are dynamically interrelated (without the artificial division of studies into the diffusion limited regime (DLR) and the surface limited regime (SLR)).

1. Introduction

Studies 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 [110]. 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 [1113], 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 [1423]. This work is a continuation of [2427], 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 [1820, 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 ( is the plate thickness and is its area). The temperature 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:where is the concentration of diffusing hydrogen. The diffusion coefficient depends on the sample temperature in an Arrhenius way with preexponential factor and activation energy .

Initial data are determined by the fact that the sample had been preliminarily degassed: , , .

Nonlinear boundary conditions are derived from the material flux balance:Here, , are the amounts of hydrogen atoms in the input chamber of volume and output chamber of volume , , . 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 is related to the pressure by the Hertz Knudsen formula: ( is the Boltzmann constant, is the mass of a hydrogen molecule). In the context of the experiment it is convenient to choose the following measurement units , . Then we numerically obtain the dependence , (, ). 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 . One can write (as a parameter of the model) instead of , but it is more convenient to interpret the dimensionless probability factor as the fraction of absorbed H atoms within the notation. Thus, 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, are the densities of the desorption flux from the sample (deviation from the square desorption law is significant only at extreme temperatures) and is the desorption coefficient. We also assume that and depend on the temperature in an Arrhenius way. Formally, the activation energy 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 is maintained at a constant temperature on both sides of the membrane, the equilibrium concentration of the dissolved atomic diffusionally mobile hydrogen is finally established. By equating the derivatives in model (2) to zero, we get Sieverts’ law , .

Let us clarify the experimental conditions. The volumes comprise several liters, the thickness of the membrane is , the area is about , and the inflow pressure is within several hundreds .

It now remains to find the magnitudes of , . Within the time of transfer through the membrane the gas is in the thermodynamical quasi-equilibrium, wherefore we use the formula . Here, is the number of gas particles occupying the volume at the temperature and the pressure (in the SI system , , ). Taking into account the relationships and (formally), we get the following expressions for the corresponding pressures and volumes in the boundary conditions (2): , . Here, , , are the numerical values in the selected units (, , ).

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 . The gas inside the volumes (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 (and the subsequent ones, only excluding the value ). The need of such a refinement arises from the characteristics of this particular experimental unit. Such an adjustment of the values of 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 .

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: Let us formulate the problem of modelling the concentrations using the pressures (the problem is also of interest per se) without the quasi-equilibrium simplification . The quasi-stationary state is achieved within a time , which is short compared to the total experiment time (). So, the original model (1)-(2) can be simplified (taking into account the formula ):Since by virtue of the “inlet-outlet” balance the equalities hold true, it is sufficient to express the concentrations from the boundary conditions (5) and substitute them into the first equation of (4) (the sign is chosen depending on whether the index is or ). The dimensionless normalized variables are convenient for numerical simulation:In addition, the system of equations (5) is compactly written in the symmetric form , . For the variable we obtain the incomplete quartic equation , 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 . Therefore, we shall aim to derive differential equations for , since information about the dynamics of the boundary concentrations is also important.

Differentiate (5) with respect to time and substitute the pressure derivatives from (4). For the variables we get the system

The change of variables in (7) determines the concentrations from which the model pressures , are calculated using (5).

Let us formulate step-by-step the numerical algorithm of modelling the pressures for the current values of , , coefficients (the authors used the Scilab freeware). We target at the “normal” experimental conditions [11, 12, 2931], including the values of , , , , .

(1) We fix : 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 we choose the root of the biquadratic polynomial . From physical considerations it follows that and thus .

(2) The system of equations , () yields the missing value of . Formally, one equation is enough, but we take into account averaging procedures including determination of the values of .

(3) We numerically integrate the ODE system (8), (9) with the obtained initial data. The change of variables in (7) defines the concentrations , which are used to calculate the model pressures , from (5).

Computational experiments show that the model curves almost coincide (at ) 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 . All the variable parameters of the original model that influence the permeability, namely, , , , 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 [1113, 2933]. We fix , , , , , . We set the value and calculate the corresponding desorption coefficient . The input pressures 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 . 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): . Within the model we determine and by the formulas For the given values of the parameters we have , , .

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 The penetrating flux is .

From the computational point of view it is convenient to introduce the dimensionless time which is oriented at the characteristic diffusion time . As small a singularity appears if we directly use the partial sum of the expressionLet us provide another expression for using the properties of the Jacobi theta function. More precisely, we are interested in the function We have an alternative presentation for [34]:The series on the left is rapidly converging for large . But the series on the right is rapidly converging for small . If we define (), then we get , or , which is known as the functional equation for the theta function [35]. After some auxiliary transformations for we get the appropriate representation

The function has an -shaped form of the saturation curve (see the inset in Figure 1). Stage I ends with .

Stage II. is the new time zero: Stage II ends with . The penetrating flux is

Stage III. The formulas are similar to the cyclic interchange , .

The result of connecting the stages into a single “experimental” curve of the penetrating flux (conventionally, , ) is shown in Figure 1.

3.2. Converting Pressure into Flux

During the real experiment, the gas pressure inside the outlet volume 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 . We take as the framework the ideal gas state equation: . Here, , , is the number of particles ( molecules), and is the Boltzmann constant. Differentiating on we get . Let us calculate the particle balance over the time : where (the concentration), . The factor is due to the fact that the diffusion flux is atomic, and the particle in the volume is the molecule. Let us divide the equation of the material balance by and direct . Finally, we get the differential equation: Denoting , , we get (), wherefrom The measured function is noised, so we first apply a smoothing procedure and only then calculate the derivative . 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 . Asymptotically, a steady state value is established: .

Only the operation of integration is needed to calculate the lag time (see below): so a preliminary approximation of the derivative is unnecessary. It suffices to use the composite Simpson formula. One should remember that we use the following measurement units within this subsection , , (), . 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 ( are measured). Considering the new time zero we also have . We are so far talking about the direct problem of modelling hydrogen pressures inside the volumes . We specify the values: , , , , . We target here the experimental conditions and the data for the alloy [11, 33].

We numerically solve the initial boundary problem of hydrogen permeability:

The membrane temperature is taken in the dependences of , , , on , and the temperature of the gas inside the chambers is taken in the expressions for (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 into the expressions for ). To this end one should skip the initial time 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 asymptotically moves to the stationary value . Hence, (). The intersection of the asymptote with the axis gives the so-called lag time , which allows estimating the diffusion coefficient. Analytically, Note that the value under the integral sign is a relative magnitude which does not require absolute values of the penetrating flux in any measurement units (). In addition, the value does not depend on . It is usually assumed that the locally equilibrium concentration is quickly established at the inlet, so one can additionally estimate the solubility and the permeability using the corresponding value . 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 () is quickly established given that . The value of as such is yet to be clarified. Thus, only the estimate of the diffusion coefficient is considered reliable at this stage.

With a new zero time reference, integrating the expression we get where , . Formally, we obtain the same expressions for the lag time and the estimate of if we change both the zero time and the flux baseline ( value overstatement). There is no additional information here (about the target values of the surface parameters and ), but the triple penetration experiment allows refining . For the model numerical experiment we have .

4.2. Isotherm: Initial Estimates of ,

In experimental practices it is common to draw and analyze the isotherm, i.e., the curve of the steady state penetrating flux dependence on the inlet pressure , while vacuum pumping is performed at the outlet side. If one targets the Sieverts’ law and the (quasi)equilibrium concentration at the inlet side (hence ), then it is natural to plot the dependence .

Let us analyze the steady state flux balance equation: , Asymptotic analysis shows that the dependence has a parabolic shape at low inlet pressures : and a straight line form at relatively high pressures : Using the straight-line segment of the isotherm we find (the slope of the straight line), and knowing the estimate of we determine the initial approximation of the solubility coefficient . From the intersection of the straight line with the ordinate axis we find . Knowing the values of and we compute and .

A graphic illustration (using a minimal required set of the calculated model 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 , , 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 , , 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 , , corresponds to this pressure . 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 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 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 . This is the bulk desorption coefficient (the coefficient of effective recombination of atoms into molecules): , .

The heating is usually linear . The heating rate is not too high (<K/s). When the temperature maximum is reached (if degassing is not yet completed), then heating is stopped: . 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: , ,Here, are the surface concentrations (); is the parameter of local equilibrium between the surface and the subsurface bulk (coefficient of quick dissolution); is the surface desorption coefficient: 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:Coefficients , are the descriptors of the rate of dissolution in the bulk and transfer to the surface. When concentrations are nowhere near maximum and on the relative scale, we obtain (29), where . If the surface is isotropic (in terms of ), then the parameter is a little dependent on temperature. The density of the H atoms adsorption flux can be modeled by the term for balance equations (30), (31). For the ranges of weak concentrations and sufficiently high working temperatures the degree of surface coverage satisfies . 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 The generalized quick dissolution coefficient of local “surface-bulk” equilibrium is assumed to be constant () in the given heating range. The TDS method fulfills the symmetry conditions:The information capacities of the initial and final stages of the TDS experiment are low. So it is sufficient to limit ourselves to , 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 or TDS spectra () for different saturation conditions and heating rates. The conversion is made more specific by taking into account the parameters of the experimental unit. Modern equipment provides a means of creating deep vacuum (). For this reason, the component () is the key control for the saturation stage, but resorption is neglected for the degassing stage ().

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 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:where is the concentration of hydrogen atoms captured by defects of different types; are the coefficients of H capture and release by traps; is the defects saturation degree . 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 is identically zero and value is positive only if the critical temperature is reached: . It is easy to simulate the required number of TDS-peaks using different binding energies (coefficients ). Numerical algorithms based on difference schemes and modeling results were presented in [24, 26]. In this paper we use only (1) (). 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 value.

(1) On the surface the following ratio is applied: Here, is the degree of surface coverage (next have a similar meaning). Let us formulate :

(2) In the “surface–bulk” equilibrium we have (, , , , then from Formula (33) we obtain

(3) For definiteness, we take into account only one type of traps (). We use (36). For symmetry, we add the saturation factor for . Here we have , ; then we obtain

The “saturation-degassing” experiment provides information about a general average concentration in the bulk (end surfaces are neglected): Normalize by and consider the dependenceLet us focus on the curves in the axes (). Numerical results are presented in Figure 5 (, , 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 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 . 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:Membranes are usually very thin and their permeability is relatively high. We can accept for the concentration gradient . In accordance with the Fick’s law, the formula is “exact”. Here, the function 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 is substituted with the equilibrium concentration in accordance with Sieverts’ law ( is the solubility coefficient, or solubility for short). On account of the inequalities , , this substitution overrates the second factor in the formula . Hence, the value of needs to be formally marked down to retain the equality. So, if the value of (permeability) is found from experimental data by fitting as indicated in formula (43), then we obtain . The and values from the (typically used) formula determine the lower bound of the diffusion coefficient . 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 , . The “saturation-degassing” experiment yields the total value and an overstated (for the permeability problem) value .

Hence, we can estimate the value by fitting using formula (43). This information has practical value as a convenient coefficient for translation from pressures to flux. If, however, the value is taken from one experiment (or paper) and the value is taken from another source, then, strictly speaking, we get a ranking of three different numbers . 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 (as a parameter of (43)) has an -form (Arrhenius-like) of the saturation curve based on the order of pressures. It is only for relatively high pressures (where are near Sieverts’ concentrations) that we get .

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 ( usually is fractions of mm). Model (36) gives the possibility to get any number of peaks using traps with different parameters . 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 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 (). It is reasonable to try to avoid iterative solution of the boundary value problem for interim approximations of the model parameters , , , , , , . 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 , , Let us replace the time (keeping the former notation ):Here 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: Let us write the solution of the linear boundary value problem using the source function (Green’s function): Boundary conditions contain the derivative : For we have divergent series. So, term-by-term integration is implied. For the original time we get , Finally, the dynamic boundary condition is written in the integrodifferential form: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 uniquely determines the solution . The analogy with the functional differential equation of the form of the neutral type [38] is that it is impossible to eliminate the derivative on the right side of the equation through integration by parts lest a divergent series arises. We are concerned with the time interval , which corresponds to the TDS peak. Measurement for 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: , , (). Retaining the notation , we obtainwhere is a dimensionless parameter of quadratic “desorption”.

6.3. Initial Data

We specify the factor for quadratic nonlinearity. Initial saturation is conducted under relatively high temperature and pressure . After the steady state of saturation is attained we get

Hereafter, let us be guided by a maximum limit of surface concentration at around (monolayer in the context of geometric statics). If during initial saturation the surface concentration for the given model parameters is , the degree of surface coverage must be taken into account. Then, to calculate the initial value of we take the ratio instead of .

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 may be higher, if it is specified what meaning is attached to the term “bubbling” surface layer (at a quite high temperature). The real 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” 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, takes finite values for . The series is rapidly converging for large . Formally, substituting (the integration variable reaches the upper limit ) we obtain a divergent series. This can be “fixed” by term-by-term integration. The order of terms in the series is (), 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): The series on the left is rapidly converging for large . But the series on the right is rapidly converging for small . Let us run auxiliary transformations: Proceeding from here, the last series is subtracted from the first series and the result is doubled. The following expression is obtained for : . The series is very rapidly converging for small . As we have and the integrable singularity . The graph for has an -shaped (“Arrhenius-like”) saturation curve form and . The function increases monotonically to a maximum (≈0.828 for ) and then decreases monotonically. To represent the series 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 we obtain the equationwhere the fraction (weak singularity under the integral) rapidly decreases from infinity () to near zero (). For , the lower limit of the integral can be replaced with . Thus, the neglectable background can be easily identified in the original physical time using the relation . The compact functional differential () TDS equation (58) with initial data 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 . The theory of integrodifferential equations (the “residual” integral containing the function 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 cm, K, . The assumed values of model parameters are the following: (steel)  cm2/s, ,  cm2/s, , , , , , at.H/cm3, K, K/s; (Ni) , , , , , , , ,  Torr, , ; (W) , , , , , , , , , , ; (Be) , , , , , , , , , , .

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 ( for dimensionless time, ). Let us fix due to the smooth function graph (Figure 4) and the trapezoid rule for numerical integration. The replacement of is naturally targeted at the diffusion time scale so that the step 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 .

The greatest contribution to the integral is made by the value , due to nonlimited (but integrable) singularity. Thus, the quadratic approximation is used due to function concavity. Let us consider the current segment of time , . The conditions determine the values of , (constants on ): . Represent the integral from Formula (58) as a sum , . For the second additive, the singularity is explicitly integrated by substituting (9). The trapezoid rule and mean value theorem are used for the integral without singularity To approximate the first integral we use only several terms (for definiteness ) on the right-hand side taking into account the factor . The negative additives () thereby drop out. Let us compensate for that by replacing with and introducing the additional variables : As a result we obtain an ODE system instead of (58):Ellipsis stands for the right-hand side of the first line of the first equation (). The series for converges very rapidly. It is sufficient to use 5–8 terms for software.

The sought function on the current segment of dimensionless time (in the system (65) , , ) 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 should be additionally used to improve the modeling accuracy. Returning to the physical time, we get the model spectrum ().

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 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 and DLR corresponds to . 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 , , , , , , . Determine the values of , , under the saturation conditions , as described in Initial saturation stage. If the equilibrium bulk concentration of dissolved diffusing hydrogen is known, then presetting of 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 length interval). Other ODE approximations were presented in [27]. This procedure is standard for every modern software.(iii)The graph or spectrum is plotted in physical time using the function .

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 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 is known (Daines-Berrer method is usually used for the DLR of permeability). Let us formulate the problem: to estimate the surface parameters , , using the (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 () for the situation where accumulation on the surface can be neglected () was presented in [25].

Since the function is known, it is reasonable to move directly to dimensionless time , which is oriented at the characteristic time of the diffusion transfer . The old notation is retained not to complicate the formulas. Let us present the nonlinear additive term in the system (65) in the following form: Here, is a function of the parameter . It is obtained through transformations using notations from (52) and (65). Elementary but somewhat lengthy formulas are omitted. The function () is known from measurements. After substitution into (65) we obtain a system of three linear ODE. It is reasonable to add the variables to improve the calculation accuracy. The measurements are usually noisy, but the function 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 . The identification algorithm can thus be decomposed. By setting the value of we can eventually calculate the dependence in the original physical time . Be reminded now thatwhere , . Since heating is monotonic , we can introduce the coordinates , . The parametric curve , is obtained on the plane . 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 : the curve , has to be a straight line segment on the plane (). This subproblem is scalar, where only varies.

Formally, we prolong the straight line segment until it crosses the coordinate axes. The crossing with the -axis () yields the value. The crossing with the abscissa determines the 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 , ) 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 yields the mass balanceThe experiments have demonstrated that the accuracy of such estimation decreases for . Curves based on the “correct” numerical value (the resultant curve is close to the straight line segment) and curves for 20% deviation from 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” value is exceeded. Mathematical singularity appears because, formally, the function 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 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 , desorption , and absorption 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 , , 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.

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).